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

(*Version April 29, 2023 of displayRREF*)
(*bug reports are welcome*)
(*updated august 20, 2023 to fix a minor bug*)

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,entry},
    displayMatrix[A,dashed];
    {nRows,nCols} = Dimensions[A];
    
    keepEliminating=True;
    
    n=1; m=1; (*n is row number, m is column number *) 
    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]];
                       A=Simplify[A];
                       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;;]];
                                  A=Simplify[A];
                                  Print@makeNiceMatrix[A,{n,m},dashed];
                            ]
                       ,{j,n+1,nRows}
                       ];
                  ];
                  pivotsFound = AppendTo[pivotsFound,{n,m}];
                  If[n==nRows,
                      keepEliminating=False
                  ,
                      n++;
                      If[m<nCols,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;
                                  A=Simplify[A];
                                  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@makeNiceMatrix[A,{n,m},dashed];
    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,;;]];
                      A=Simplify[A];
                      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}}]}}];

mat={{1,2},{3,4},{5,6}}
{result,pivots}=displayRREF[mat]

(*example usages*)
mat={{1,2,3},{4,5,6}}
{result,pivots}=displayRREF[mat]