> with(linalg);
Consider the following augmented matrix of a system of equations. We need to transform it into the reduced row-echelon form by the row operations.

> A:=matrix(5,8,[[4,8,-2,-5,0,-23,2,37],[1,2,2,5,-5,13,-2,3],[2,4,-2,-5,2,-19,2,21],[3,6,1,0,0,-1,0,22],[4,8,1,2,-5,-1,-1,30]]);

A := [ 4 8 -2 -5 0 -23 2 37 ]
[ 1 2 2 5 -5 13 -2 3 ]
[ 2 4 -2 -5 2 -19 2 21 ]
[ 3 6 1 0 0 -1 0 22 ]
[ 4 8 1 2 -5 -1 -1 30 ]

> B:=mulrow(A,1,1/4);

B := [ 1 2 -1/2 -5/4 0 -23/4 1/2 37/4 ]
[ 1 2 2 5 -5 13 -2 3 ]
[ 2 4 -2 -5 2 -19 2 21 ]
[ 3 6 1 0 0 -1 0 22 ]
[ 4 8 1 2 -5 -1 -1 30 ]

> B:=addrow(B,1,2,-1);

B := [ 1 2 -1/2 -5/4 0 -23/4 1/2 37/4 ]
[ 0 0 5/2 25/4 -5 75/4 -5/2 -25/4 ]
[ 2 4 -2 -5 2 -19 2 21 ]
[ 3 6 1 0 0 -1 0 22 ]
[ 4 8 1 2 -5 -1 -1 30 ]

> B:=addrow(B,1,3,-2);

B := [ 1 2 -1/2 -5/4 0 -23/4 1/2 37/4 ]
[ 0 0 5/2 25/4 -5 75/4 -5/2 -25/4 ]
[ 0 0 -1 -5/2 2 -15/2 1 5/2 ]
[ 3 6 1 0 0 -1 0 22 ]
[ 4 8 1 2 -5 -1 -1 30 ]

> B:=addrow(B,1,4,-3);

B := [ 1 2 -1/2 -5/4 0 -23/4 1/2 37/4 ]
[ 0 0 5/2 25/4 -5 75/4 -5/2 -25/4 ]
[ 0 0 -1 -5/2 2 -15/2 1 5/2 ]
[ 0 0 5/2 15/4 0 65/4 -3/2 -23/4 ]
[ 4 8 1 2 -5 -1 -1 30 ]

> B:=addrow(B,1,5,-4);

B := [ 1 2 -1/2 -5/4 0 -23/4 1/2 37/4 ]
[ 0 0 5/2 25/4 -5 75/4 -5/2 -25/4 ]
[ 0 0 -1 -5/2 2 -15/2 1 5/2 ]
[ 0 0 5/2 15/4 0 65/4 -3/2 -23/4 ]
[ 0 0 3 7 -5 22 -3 -7 ]

> B:=mulrow(B,2,2/5);

B := [ 1 2 -1/2 -5/4 0 -23/4 1/2 37/4 ]
[ 0 0 1 5/2 -2 15/2 -1 -5/2 ]
[ 0 0 -1 -5/2 2 -15/2 1 5/2 ]
[ 0 0 5/2 15/4 0 65/4 -3/2 -23/4 ]
[ 0 0 3 7 -5 22 -3 -7 ]

> B:=addrow(B,2,3);

B := [ 1 2 -1/2 -5/4 0 -23/4 1/2 37/4 ]
[ 0 0 1 5/2 -2 15/2 -1 -5/2 ]
[ 0 0 0 0 0 0 0 0 ]
[ 0 0 5/2 15/4 0 65/4 -3/2 -23/4 ]
[ 0 0 3 7 -5 22 -3 -7 ]

> B:=addrow(B,2,4,-5/2);

B := [ 1 2 -1/2 -5/4 0 -23/4 1/2 37/4 ]
[ 0 0 1 5/2 -2 15/2 -1 -5/2 ]
[ 0 0 0 0 0 0 0 0 ]
[ 0 0 0 -5/2 5 -5/2 1 1/2 ]
[ 0 0 3 7 -5 22 -3 -7 ]

> B:=addrow(B,2,5,-3);

B := [ 1 2 -1/2 -5/4 0 -23/4 1/2 37/4 ]
[ 0 0 1 5/2 -2 15/2 -1 -5/2 ]
[ 0 0 0 0 0 0 0 0 ]
[ 0 0 0 -5/2 5 -5/2 1 1/2 ]
[ 0 0 0 -1/2 1 -1/2 0 1/2 ]

> B:=swaprow(B,3,4);

B := [ 1 2 -1/2 -5/4 0 -23/4 1/2 37/4 ]
[ 0 0 1 5/2 -2 15/2 -1 -5/2 ]
[ 0 0 0 -5/2 5 -5/2 1 1/2 ]
[ 0 0 0 0 0 0 0 0 ]
[ 0 0 0 -1/2 1 -1/2 0 1/2 ]

> B:=mulrow(B,3,-2/5);

B := [ 1 2 -1/2 -5/4 0 -23/4 1/2 37/4 ]
[ 0 0 1 5/2 -2 15/2 -1 -5/2 ]
[ 0 0 0 1 -2 1 -2/5 -1/5 ]
[ 0 0 0 0 0 0 0 0 ]
[ 0 0 0 -1/2 1 -1/2 0 1/2 ]

> B:=addrow(B,3,5,1/2);

B := [ 1 2 -1/2 -5/4 0 -23/4 1/2 37/4 ]
[ 0 0 1 5/2 -2 15/2 -1 -5/2 ]
[ 0 0 0 1 -2 1 -2/5 -1/5 ]
[ 0 0 0 0 0 0 0 0 ]
[ 0 0 0 0 0 0 -1/5 2/5 ]

> B:=swaprow(B,4,5);

B := [ 1 2 -1/2 -5/4 0 -23/4 1/2 37/4 ]
[ 0 0 1 5/2 -2 15/2 -1 -5/2 ]
[ 0 0 0  1 -2 1 -2/5 -1/5 ]
[ 0 0 0 0 0 0 -1/5 2/5 ]
[ 0 0 0 0 0 0 0 0 ]

> B:=mulrow(B,4,-5);

B := [ 1 2 -1/2 -5/4 0 -23/4 1/2 37/4 ]
[ 0 0 1 5/2 -2 15/2 -1 -5/2 ]
[ 0 0 0 1 -2 1 -2/5 -1/5 ]
[ 0 0 0 0 0 0 1 -2 ]
[ 0 0 0 0 0 0 0 0 ]
Now the matrix is in the row echelon form but not yet in the reduced row echelon form. So we have to continue.

> B:=addrow(B,2,1,1/2);

B := [ 1 2 0 0 -1 -2 0 0 ]
[ 0 0 1 5/2 -2 15/2 -1 -5/2 ]
[ 0 0 0 1 -2 1 -2/5 -1/5 ]
[ 0 0 0 0 0 0 1 -2 ]
[ 0 0 0 0 0 0 0 0 ]

> B:=addrow(B,3,2,-5/2);

B := [ 1 2 0 0 -1 -2 0 0 ]
[ 0 0 1 0 3 5 0 -2 ]
[ 0 0 0 1 -2 1 -2/5 -1/5 ]
[ 0 0 0 0 0 0 1 -2 ]
[ 0 0 0 0 0 0 0 0 ]

> B:=addrow(B,4,3,2/5);

B := [ 1 2 0 0 -1 -2 0 0 ]
[ 0 0 1 0 3 5 0 -2 ]
[ 0 0 0 1 -2 1 0 -1 ]
[ 0 0 0 0 0 0 1 -2 ]
[ 0 0 0 0 0 0 0 0 ]
That's it. This matrix is in the row echelon form. The leading unknowns are x1, x3, x4, x7, the free unknowns are x2, x5, x6. The general solution is: 
                              x1=8-2s+t+2u 
                              x3=-2-3t-5u 
                              x4=-1+2t-u 
                              x7=-2 
                              x2=s 
                              x5=t 
                              x6=u