We first express the system of equations as a single matrix equation,

(1)Mx=b

Since there are as many equations (rows in M) as there are variables (columns in M), M will be a square matrix. If there is a solution to this set of equations, then there will exist another matrix M1 such that M1M=I, where I is the identity matrix. In that case,

(2)M1Mx=M1b
(3)Ix=x=M1b

There are algorithmic methods for finding M1 (and thus x; in our case we’ll be using standard libraries to calculate this matrix and the corresponding solution. One standard method is to have the computer calculate M1 directly using an “invert()” function, and then multiply that inverse matrix by b to obtain x.

Most computational packages that have the ability to calculate M1 also have an additional function, often called “solve()”, that takes care of the multiplication as well, and gives the user x directly.

For more information about how these functions work, see Numerical Recipes in C by Flannery, Teukolsky, Press, and Vetterling, or similar text.