(*any bugs, please report to Nasser M. Abbasi  nma@12000.org  so I can correct the code *)

SetDirectory[NotebookDirectory[]]
CurrentValue[$FrontEnd,{"PrintAction"}]={"PrintToNotebook"};
CurrentValue[$FrontEnd,{"MessageOptions","KernelMessageAction"}]={"PrintToNotebook","Beep"};

displayMatrix[A_?(MatrixQ[#]&),dashed_?BooleanQ]:= If[dashed,displayMatrix[A,Last@Dimensions[A]],Print[MatrixForm[A]]]

displayMatrix[A_?(MatrixQ[#]&),nCols_Integer]:= Print[A[[All,1;;nCols]]//matWithDiv[nCols,Background->LightOrange]]

displayRREF[Ain_?(MatrixQ[#]&),dashed_:False]:=Module[
  {multiplier,j,i,pivotRow,pivotCol,nRows,nCols,p,tmp,startAgain,A=Ain,n,m,pivotsFound={},keepEliminating,nIter},
    displayMatrix[A,dashed];
    {nRows,nCols} = Dimensions[A];
    
    keepEliminating=True;
    
    n=1; m=1;
    nIter = 0;
    While[keepEliminating,
        nIter++;
        If[nIter>100, (*safe guard*)
             Return["Internal error. Something went wrong. Or very large system?",Module]
        ];
        
        If[dashed && m==nCols,
             keepEliminating=False 
        ,
             Print["Pivot is A(",n,",",m,")"];
             Print@makeNiceMatrix[A,{n,m},dashed];
             If[A[[n,m]] != 0,
                  If[A[[n,m]] != 1,
                       A[[n,All]] = A[[n,All]]/A[[n,m]];
                       Print["Making the pivot 1 using using row(",n,")= row(",n,")/A(",n,",",m,")"]
                       (*Print@makeNiceMatrix[A,{n,m},dashed]*)
                 ];
                 If[n<nRows,
                       Do[
                            If[A[[j,m]] != 0,
                                  multiplier = A[[j,m]]/A[[n,m]];
                                  Print["Zeroing out element A(",j,",",m,") using row(",j,")=",multiplier,"*row(",n,")-row(",j,")"];
                                  A[[j,m;;]] = A[[j,m;;]] - multiplier*A[[n,m;;]];
                                  Print@makeNiceMatrix[A,{n,m},dashed];
                            ]
                       ,{j,n+1,nRows}
                       ];
                  ];
                  pivotsFound = AppendTo[pivotsFound,{n,m}];
                  If[n==nRows,
                         keepEliminating=False
                  ,
                      n++;
                      m++   
                  ]
            ,
                  (*pivot is zero*)
                  Print["Pivot is zero"];
                  (*Print@makeNiceMatrix[A,{n,m},dashed];*)
                  If[n==nRows&&m==nCols, keepEliminating=False
                  ,
                       (*pivot is zero. If we can find non-zero pivot row below, then exchange rows*)
                       If[n<nRows,
                             p = FirstPosition[A[[n+1;;,m]],_?(#!=0&)];
                             If[ p===Missing["NotFound"]|| Length[p]==0,
                                  If[m<nCols,
                                       m++
                                  ,
                                       keepEliminating=False
                                  ]
                             ,
                                  (*found non zero pivot below. Exchange rows*)
                                  tmp = A[[n,All]];
                                  A[[n,All]] = A[[First[p]+n,All]];
                                  A[[First[p]+n,All]] = tmp;
                                  Print["Exchanging row(",n,") and row(",First[p]+n,")"]
                                  (*Print@makeNiceMatrix[A,{n,m},dashed]*)
                            ]
                      ,
                            If[m<nCols,
                                  m++
                            ,
                                  keepEliminating=False
                            ]
                      ]
                  ]
             ]
        ]
    ];
         
    
    (*pivotsFound = DeleteDuplicates[pivotsFound];*)
    Print[">>>>>>Starting backward elimination phase. The pivots are ",pivotsFound];
    Do[ 
        pivotRow=First@entry;
        pivotCol=Last@entry;
        
       
        If[pivotRow>1,
            Do[
                 If[ A[[i,pivotCol]] != 0,
                      Print["Zeroing out element A(",i,",",pivotCol,") using row(",i,")=row(",i,")-A(",i,",",pivotCol,")*row(",pivotRow,")"];
                      A[[i,;;]] = A[[i,;;]] - A[[i,pivotCol]]*A[[pivotRow,;;]];
                      Print@makeNiceMatrix[A,{pivotRow,pivotCol},dashed]
                 ]
            ,
            {i,pivotRow-1,1,-1} 
           ]
        ]
        ,
        {entry,pivotsFound}    
   ];
   
   {A,pivotsFound[[All,2]]}
]

makeSolutionSpecialCase[A_?(MatrixQ[#]&),b_?(VectorQ[#]&),pivotCols_List]:=Module[
   {nRows,nCols,nLeadingVariables,nFreeVariables,n,m,k,variables={},eq,freeVariables,sol={}},
   
    Print["Pivot columns are ",MatrixForm[pivotCols]];
    
    ClearAll[x,t]; (*did not make them local, to prevent $ from showing in print*)
    {nRows,nCols} = Dimensions[A];
    nLeadingVariables = Length[pivotCols];
    nFreeVariables =  nCols-nLeadingVariables;    
    Print["There are ",nLeadingVariables," leading variables and ",nFreeVariables," free variables. These are "];
    
    Array[t, nFreeVariables];
    Array[x, nCols];
    
    m=0;k=0;
    Do[
        If[Not[MemberQ[pivotCols,n]],
            m++;
            Print[x[n]," is a free variable. Let ",x[n],"=",t[m]];
            AppendTo[variables,t[m]];
            AppendTo[sol,0]
        ,
            Print[x[n]," is a leading variable"];
            AppendTo[variables,x[n]];
            AppendTo[sol,b[[++k]]]
        ]
    ,{n,1,nCols}
    ];
    
    freeVariables=(t[#]&/@Range[nFreeVariables]);
    Print["Hence the system after RREF is the following>>>>>>"];
    Print[MatrixForm[A.variables],"=",MatrixForm[b]];
    Print["There is different solution for different value of the free variables."];
    Print["Setting free variable ", freeVariables," to zero gives"];
    variables=variables/.((t[#]->0)&/@Range[nFreeVariables]);
    Print[MatrixForm[A.variables],"=",MatrixForm[b]];
    Print["Therefore the final solution is "];
    Print[MatrixForm[x[#]&/@Range[nCols]],"=",MatrixForm[sol]]
    
   
]

(*version 3/10/2017. Original version*)
(*version 5/7/2022. Make it handle all cases*)
(*version 5/12/2022. Added frame around pivot as it moves*)

displayRREF[A_?(MatrixQ[#]&),b_?(VectorQ[#]&)]:=Module[
    {B,(*augmented*)nRows,nCols,nEquations,nVariables = Length@b,rref,pivotCols,matRank,augmentedRank},
    
    {nRows, nCols} = Dimensions[A];
    nEquations     = nRows;
   
    If[nEquations != nVariables, Return["Size of b vector is not the same as number of rows in the matrix",Module]];
   
    matRank   = MatrixRank[A];
    B = Join[A,Transpose[{b}],2];
    
    Print["Augmented matrix is "];
    displayMatrix[B,True];

    augmentedRank = MatrixRank[B];

    If[matRank<augmentedRank, (*Case No solution*)
         Print["System is not consistent, no solution exist. Try using least squares."]
    ,
        (*we must have rank A == rank[A|b]  -- system is consistent. It can have*)
        (*infinite solutions or one unique solution*)
        If[matRank == nCols, 
             If[nCols == nRows, (*square. case A*)
                  Print["System is consistent. One unique solution"];
                  {rref,pivotCols} = displayRREF[B,True];
                  
                  If[MatrixRank[A]!=Length[pivotCols], (*Verify*)
                      Print["Internal Error detected. Pivot columns not same as Rank. Please report bug"]
                  ,
                      Print["Solution  vector is ", MatrixForm[rref[[All, nCols + 1]]]]
                  ]
            ,
                Print["Internal Error detected. matRank != nCols but not square. "]
            ]
       ,
            (*-- rank A  < N. Case C *)
            Print["System is consistent but infinite solutions."];
            {rref,pivotCols} = displayRREF[B,True];
            makeSolutionSpecialCase[rref[[All,1;;nCols]],rref[[All,nCols+1]],pivotCols]
      ] 
   ]
]

makeNiceMatrix[mat_?MatrixQ,pivot_List,dashed_?BooleanQ]:=Module[{g,nRow,nCol},
    {nRow,nCol} = Dimensions[mat];
   (*g=Grid[mat,Dividers->{dashPosition->{Red,Dashed}},Background->{None,None,{pivot->Pink}}];*)
   If[dashed,
       g = Grid[mat,Dividers->{nCol->{Red,Dashed}},Frame->{None,None,{pivot->True}},Background->LightOrange]
   ,
       g = Grid[mat,Frame->{None,None,{pivot->True}}]
   ];
   MatrixForm[{{g}}]
]

(*thanks to http://mathematica.stackexchange.com/questions/60613/how-to-add-a-vertical-line-to-a-matrix*)
(*makes a dash line inside Matrix*)
Format[matWithDiv[n_,opts:OptionsPattern[Grid]][m_?MatrixQ]]:=MatrixForm[{{Grid[m,opts,Dividers->{n->{Red,Dashed}}]}}];

Example 0

A={{1,-2,4,3,2},{-7,14,-28,-12,-23},{4,-8,12,12,8},{-3,6,-4,-1,-14}};
b={2,4,18,-10};
displayRREF[A,b]
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]

Example 1

A= {{2,-2,4,3,4},{-7,14,-28,-12,5},{4,-8,12,12,6}};
b={2,4,18};
displayRREF[A,b]
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]

Example 2

(*no solution example*)
A= {{1,-2,-1,3},{-2,4,5,-5},{3,-6,-6,8}};
b={1,3,1};
displayRREF[A,b];
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]

Example 3

A= {{2,1,4},{2,-3,4},{3,-2,6}};
{rref,p}=displayRREF[A];
Print["pivots columns are",MatrixForm[p]]

Example 4

A= {{1,2,1,1},{3,5,1,3},{2,6,7,1}};
{rref,p}=displayRREF[A];
Print["pivots columns are",MatrixForm[p]]

Example 5

A= {{1,-2,-1,3,1},{-2,4,5,-5,3},{3,-6,-6,8,1}};
{rref,p}=displayRREF[A];
Print["pivots columns are",MatrixForm[p]]

Example 6

A= {{6,20},{2,7}};
b={-8,4};
displayRREF[A,b]
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]

Example 7

A= {{-6,4},{3,-2}};
b={2,-1};
displayRREF[A,b]
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]

Example8

A= {{2,5},{10,25}};
b={10,50};
displayRREF[A,b]
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]

Example 9

A= {{-1,1,0},{1,-2,1},{0,1,-1}};
b={0,0,0};
displayRREF[A,b]
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]

Example 10

A= {{2,-1,2,-1,3},{1,2,3,1,0},{1,0,-2,0,-5}};
b={14,15,-10};
displayRREF[A,b]
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]

Example 11

A= {{1,3,-1},{4,-1,2},{2,-1,-3}};
b={-4,3,1};
displayRREF[A,b]
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]

Example 12

A= {{1,3,-1},{4,-1,2},{2,-7,4}};
b={4,8,0};
displayRREF[A,b]
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]

Example 13

A= {{-3,-5,36},{-1,0,7},{1,1,-10}};
b={10,5,-4};
displayRREF[A,b]
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]

Example 14

A= {{1,1,1},{2,3,2},{-2,-1,-2}};
b={0,-3,1};
displayRREF[A,b]
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]

Example 15

A= {{0,-1,1},{0,1,0},{0,0,1}};
b={0,0,0};
displayRREF[A,b]
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]

Example 15

A= {{0,0,0},{0,0,0},{0,0,1}};
b={0,0,0};
displayRREF[A,b]
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]

A= {{-3,4,0},{5,0,9},{0,1,6}};
displayRREF[A]

A={{1,-2,4,3,2},{-7,14,-28,-12,-23},{4,-8,12,12,8},{-3,6,-4,-1,-14}};b={2,4,18,-10};
displayRREF[A,b]
Print["Mathematica gives the solution as "];MatrixForm[LinearSolve[A,b]]