Or, if you don't have a flashy calculator...

Given a system of linear equations:

a1,1x1 + a1,2x2 ... + a1,n = b1
a2,1 + ... = b2
am,1 + ... = bm

We can represent this system as a matrix equation:

A . X = B

Where A is the m*n coefficient matrix, whose (i,j) entry is ai,j,
X is [ x1, x2, ... , xn ] ,
and B is [ b1, b2, ... , bm ]

Then, we define the augmented matrix of this system to be [ A|B ]
(i.e. an m*(n+1) matrix given by adding B as an extra column to the right side of A)

This matrix can now be transformed, using elementary row operations, into row reduced echelon form.

Given a row-reduced matrix E, we can then use the Gauss-Jordan procedure to find solutions.

Case 1: The last non-zero row of E is 0,0,...,1. In this case, the system is inconsistent.

Case 2: E has n non-zero rows. Then (x1 , x2, ... , xn) = the rightmost column of E.

Case 3: E has k non-zero rows. Assign parametric values to the n-k variables whose rows in E are zero.