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

Since there are as many equations (rows in M$\mathbf{M}$) as there are variables (columns in M$\mathbf{M}$), M$\mathbf{M}$ will be a square matrix. If there is a solution to this set of equations, then there will exist another matrix M−1${\mathbf{M}}^{-1}$ such that M−1M=I${\mathbf{M}}^{-1}\mathbf{M}=\mathbf{I}$, where I$\mathbf{I}$ is the identity matrix. In that case,

There are algorithmic methods for finding M−1${\mathbf{M}}^{-1}$ (and thus x$\mathbf{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 M−1${\mathbf{M}}^{-1}$ directly using an “invert()” function, and then multiply that inverse matrix by b$\mathbf{b}$ to obtain x$\mathbf{x}$.

Most computational packages that have the ability to calculate M−1${\mathbf{M}}^{-1}$ also have an additional function, often called “solve()”, that takes care of the multiplication as well, and gives the user x$\mathbf{x}$ directly.

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