diff --git a/Driver.f90 b/Driver.f90
index fd93cbe..b0460d8 100644
--- a/Driver.f90
+++ b/Driver.f90
@@ -35,9 +35,11 @@ program Driver
 !-------------------------------------!
   integer, parameter      :: n2 = 3
   real, dimension(n2, n2) :: a2, l2, u2, p2
+  real, dimension(n2)     :: b2, x2, y2
   data (a2(1,col), col=1,n2) / 60.0, 30.0, 20.0 /  ! =--> row 1
   data (a2(2,col), col=1,n2) / 30.0, 20.0, 15.0 /  ! =--> row 2
   data (a2(3,col), col=1,n2) / 20.0, 15.0, 12.0 /  ! =--> row 3
+  data (b2(row),   row=1,n2) /  5.0, 16.0, 15.0 /
 
 !-----------------------------------+
 ! Matrix 3 for Forward substitution !
@@ -75,27 +77,22 @@ program Driver
   !   Demonstrate Gussian elimination   !
   !-------------------------------------!
   if(choice == 1) then
-    write(*,*) "Original matrix a1:"
-    call Print_Matrix(a1)
 
-    call Gaussian_Elimination(g1, b1, a1)
-
-    write(*,*) "Matrix g1 after elimination:"
-    call Print_Matrix(a1)
-    call Print_Matrix(g1)
+    ! Just print original matrix
+    call Print_Matrix("Original matrix a1:", a1)
 
-    write(*,*) "Source term vector b1 after elimination:"
-    call Print_Vector(b1)
+    ! Perform gauissian elimination on matrix and r.h.s. vector
+    call Gaussian_Elimination(g1, b1, a1)
+    call Print_Matrix("Matrix g1 after elimination:", g1)
+    call Print_Vector("R.h.s vector b1 after elimination:", b1)
 
+    ! Perform backward substitution
     call Backward_Substitution(x1, g1, b1)
+    call Print_Vector("Solution vector x1: ", x1)
 
-    write(*,*) "Solution vector x1:"
-    call Print_Vector(x1)
-
+    ! Multiply original matrix with solution vector to check result
     call Matrix_Vector_Multiply(y1, a1, x1)
-
-    write(*,*) "Vector y1 should recover the source term vector:"
-    call Print_Vector(y1)
+    call Print_Vector("Vector y1 should recover the source term:", y1)
   end if
 
   !----------------------------------------!
@@ -103,44 +100,51 @@ program Driver
   !----------------------------------------!
   if(choice == 2) then
   l2 = 0
-    call Cholesky_Factorization(l2, a2)
 
-    write(*,*) "Matrix l2 after Cholesky factorization:"
-    call Print_Matrix(l2)
+    ! Just print original matrix
+    call Print_Matrix("Original matrix a1:", a2)
 
-    call Transpose_Matrix(u2, l2)
+    ! Perform Cholesky factorization on the matrix to fin the lower one
+    call Cholesky_Factorization(l2, a2)
+    call Print_Matrix("Matrix l2 after Cholesky factorization", l2)
 
-    write(*,*) "Matrix u2 after inversion:"
-    call Print_Matrix(u2)
+    ! Transpose the lower matrix to get the upper
+    call Transpose_Matrix(u2, l2)
+    call Print_Matrix("Matrix u2 after transpose:", u2)
 
+    ! Multiply lower and upper to check
     call Matrix_Matrix_Multiply(p2, l2, u2)
+    call Print_Matrix("Matrix p2 after multiplication of l2 and u2:", p2)
+
+    ! Compute y by forward substitution
+    call Forward_Substitution(y2, l2, b2)
+    call Print_Vector("Vector y2 after forward substitution:", y2) 
 
-    write(*,*) "Matrix p2 after multiplication:"
-    call Print_Matrix(p2)
+    ! Compute x by backward substitution
+    call Backward_Substitution(x2, u2, y2)
+    call Print_Vector("Vector x2 after forward substitution:", x2) 
+
+    ! Multiply original matrix with solution vector to check result
+    call Matrix_Vector_Multiply(y2, a2, x2)
+    call Print_Vector("Vector y2 should recover the source term:", y2)
   end if
 
   !---------------------------------------!
   !   Demonstrate Backward Substitution   !
   !---------------------------------------!
   if(choice == 3) then
-    write(*,*) "Matrix a3"                          
-    call Print_Matrix(a3)
-
+    call Print_Matrix("Matrix a3:", a3)
     call Forward_Substitution(y3, a3, b3)
-
-    call Print_Vector(y3) 
+    call Print_Vector("Vector y3", y3) 
   end if
  
   !--------------------------------------!
   !   Demonstrate Forward Substitution   !
   !--------------------------------------!
   if(choice == 4) then
-    write(*,*) "Matrix a4"                          
-    call Print_Matrix(a4)
-
+    call Print_Matrix("Matrix a4:", a4)
     call Backward_Substitution(y4, a4, b4)
-
-    call Print_Vector(y4) 
+    call Print_Vector("Vector y4:", y4) 
   end if
  
   end program Driver
diff --git a/Print_Matrix.f90 b/Print_Matrix.f90
index 93cb688..c74f363 100644
--- a/Print_Matrix.f90
+++ b/Print_Matrix.f90
@@ -1,16 +1,19 @@
 !======================================================================!
-  subroutine Print_Matrix(a)
+  subroutine Print_Matrix(message, matrix)
 !----------------------------------------------------------------------!
   implicit none
 !----------------------------------------------------------------------!
-  real, dimension(:,:) :: a
+  character(len=*)          :: message
+  real,      dimension(:,:) :: matrix
 !----------------------------------------------------------------------!
   integer :: row, col  ! row used to be "i", col used to be "j"
 !----------------------------------------------------------------------!
 
-  do row = 1, size(a, 1)
-    do col = 1, size(a, 2)
-      write(*,"(f8.3)",advance="no") a(row,col)
+  write(*,*) message
+
+  do row = 1, size(matrix, 1)
+    do col = 1, size(matrix, 2)
+      write(*,"(f8.3)",advance="no") matrix(row,col)
     end do       
     write(*,*) ""
   end do           
diff --git a/Print_Matrix.int b/Print_Matrix.int
index 4b07017..97d7b67 100644
--- a/Print_Matrix.int
+++ b/Print_Matrix.int
@@ -1,6 +1,7 @@
 interface
-  subroutine Print_Matrix(a)
+  subroutine Print_Matrix(message, matrix)
     implicit none
-    real, dimension(:,:) :: a
+    character(len=*)     :: message
+    real, dimension(:,:) :: matrix
   end subroutine Print_Matrix
 end interface
diff --git a/Print_Vector.f90 b/Print_Vector.f90
index 0119149..52edd52 100644
--- a/Print_Vector.f90
+++ b/Print_Vector.f90
@@ -1,15 +1,18 @@
 !======================================================================!
-  subroutine Print_Vector(a)
+  subroutine Print_Vector(message, vector)
 !----------------------------------------------------------------------!
   implicit none
 !----------------------------------------------------------------------!
-  real, dimension(:) :: a
+  character(len=*)   :: message
+  real, dimension(:) :: vector
 !----------------------------------------------------------------------!
-  integer :: row, col
+  integer :: row
 !----------------------------------------------------------------------!
 
-  do row = 1, size(a)
-    write(*,"(f8.3)") a(row)
+  write(*,*) message
+
+  do row = 1, size(vector)
+    write(*,"(f8.3)") vector(row)
   end do           
 
   end subroutine Print_Vector
diff --git a/Print_Vector.int b/Print_Vector.int
index c070d0a..b5cff19 100644
--- a/Print_Vector.int
+++ b/Print_Vector.int
@@ -1,6 +1,7 @@
 interface
-  subroutine Print_Vector(a)
+  subroutine Print_Vector(message, vector)
     implicit none
-    real, dimension(:) :: a
+    character(len=*)   :: message
+    real, dimension(:) :: vector
   end subroutine Print_Vector
 end interface