Gaussian elimination/back-solving "by hand"

Example 1

a = {{1, 5, 2, 3}, {3, 2, 3, 5}, {-4, 3, 2, 1}}

a // MatrixForm

a[[2]] = a[[2]] - 3  a[[1]]

a // MatrixForm

a[[3]] = a[[3]] + 4  a[[1]]

a // MatrixForm

a[[3]] = 
      a[[3]] - a[[3, 2]]/a[[2, 2]]*a[[2]]

a // MatrixForm

x3 = a[[3, 4]]/a[[3, 3]]

x2 = (a[[2, 4]] - a[[2, 3]]*x3)/a[[2, 2]]

x1 = (a[[1, 4]] - a[[1, 3]]*x3 - a[[1, 2]]*x2)/
        a[[1, 1]]

RowReduce[a]

Example 2

a = {{1, 5, 2, 3}, {3, 15, 3, 5}, {-4, 3, 2, 1}}

a // MatrixForm

a[[2]] = a[[2]] - 3  a[[1]]

a // MatrixForm

a[[3]] = a[[3]] + 4  a[[1]]

a // MatrixForm

temprow = a[[2]]

a[[2]] = a[[3]]

a[[3]] = temprow

a // MatrixForm

x3 = a[[3, 4]]/a[[3, 3]]

x2 = (a[[2, 4]] - a[[2, 3]]*x3)/a[[2, 2]]

x1 = (a[[1, 4]] - a[[1, 3]]*x3 - a[[1, 2]]*x2)/
        a[[1, 1]]

RowReduce[a]

Example 3

a = {{1, 5, 2, 3}, {3, 15.0000003, 3, 5}, {-4, 3, 2, 1}}

a // MatrixForm

a[[2]] = a[[2]] - 3  a[[1]]

a // MatrixForm

a[[3]] = a[[3]] + 4  a[[1]]

a // MatrixForm

a[[3]] = 
      a[[3]] - a[[3, 2]]/a[[2, 2]]*a[[2]]

a // MatrixForm

x3 = a[[3, 4]]/a[[3, 3]]

x2 = (a[[2, 4]] - a[[2, 3]]*x3)/a[[2, 2]]

x1 = (a[[1, 4]] - a[[1, 3]]*x3 - a[[1, 2]]*x2)/
        a[[1, 1]]

RowReduce[a]

Automated Gaussian elimination/back-solving

a = {{1, 5, 2, 3}, {3, 2, 3, 5}, {-4, 3, 2, 1}}; 
Print[a//MatrixForm];
n = Length[a]; 
(* Gaussian elimination *)
Do[
    (* Check if pivot if needed *)
    p = i; 
	While[a[[p,i]] == 0, p = p + 1]; 
	If[p==n+1,Return[Print["unique solution not guaranteed"]]];
	(* Pivot if needed *)
	If[p != i, 
	   	temprow = a[[p]]; 
	    a[[p]] = a[[i]]; 
        a[[i]] = temprow];
    (* After any pivot, do elimination *) 
    Do[m = a[[j,i]]/a[[i,i]]; 
    a[[j]] = a[[j]] - m*a[[i]], {j, i + 1, n}];
    Print[a//MatrixForm],
  {i, 1, n - 1}];
(* Back solving *)
If[a[[n,n]] == 0, Return["singular matrix"]]; 
x = {}; 
x = Append[x, a[[n,n + 1]]/a[[n,n]]]; 
Do[x = Prepend[x, (a[[i,n + 1]] - Sum[a[[i,j]]*x[[j-n-1]], 
       {j, i + 1, n}])/a[[i,i]]], {i, n - 1, 1, -1}];
x

RowReduce[a]