You write a program that does Gaussian elimination. The article starts with three equations,

Code:

2x + y - z = 8
-3x - y + 2z = -11
-2x + y + 2z = -3

The augmented matrix is

Code:

2 1 -1 8 (L1)
-3 -1 2 -11 (L2)
-2 1 2 -3 (L3)

I've added the row vector name on the right, so you can easily follow what is done to each row vector, and how you get the row echelon rows.

The idea is to go over each column, from left to right. Ignoring rows already in row echelon form, select the row with the largest value in magnitude in that column. Add the row vector to all the other row vectors not yet in row echelon form, scaling so that the value in that column will be zero. Divide the selected row by its value in that column, and that row is then the next row in row echelon form.

In this case, the largest value in magnitude in the leftmost column is on the second line. That will become the first row in row echelon form (after dividing by -3, since it has that value in the first column). First, however, we'll need to add that row vector (scaled by 2/3) to the first row, and (scaled by -2/3) to the third row:

Code:

1 1/3 -2/3 11/3 (row echelon form, row R1 = -1/3 L2)
0 1/3 1/3 2/3 (L4, L4 = L1 + 2/3 L2)
0 5/3 2/3 13/3 (L5, L5 = L3 - 2/3 L3)

For the second column, the largest value in magnitude is in the last row. When scaled by 3/4, it will become the second row in row echelon form. Before that, though, we need to add the row (scaled by -(1/3)/(5/3) = -1/5) to the last row left:

Code:

1 1/3 -2/3 11/3 (row echelon form, row R1)
0 1 2/5 13/5 (row echelon form, row R2 = 3/5 L5)
0 0 1/5 -1/5 (L6, L6 = L4 - 1/5 L5)

The last row obviously has the largest value in magnitude (1/5, so scale by 1/(1/5) = 5 to get the row echelon row vector). We get

Code:

1 1/3 -2/3 11/3 (row echelon form, row R1)
0 1 2/5 13/5 (row echelon form, row R2)
0 0 1 -1 (row echelon form, row R3 = 5 L6)

The augmented matrix is now in row echelon form. It has also reduced the equations into

Code:

x + 1/3 y - 2/3 z = 11/3
y + 2/5 z = 13/5
z = -1

You can do a pass over the rows, starting at the final row, to use the augmented matrix to calculate the value for each variable:

Code:

z = -1
y + 2/5 z = 13/5
=> y = 13/5 - 2/5 z
=> y = 13/5 + 2/5 = 15/5
=> y = 3
x + 1/3 y - 2/3 z = 11/3
=> x = 11/3 - 1/3 y + 2/3 z
=> x = 11/3 - 1 - 2/3
=> x = 2