Let c be a column vector and solve M*c = [1,0,0,0,…0]. This is a system of simultaneous equations that can be solved via gaussian elimination, as described earlier. Since M is nonsingular, the solution is unique. When you're done, M*c gives the first column of the identity matrix.
Solve another system of linear equations, M*c = [0,1,0,0,0,…0]. This is the second column of R, i.e. M*c gives the second column of the identity matrix. Do this for every column in the identity matrix and build R.
Actually there is a way to run gaussian elimination once and derive all the columns of R in parallel. Place a copy of the identity matrix to the right of M. Now our augmented matrix has n rows and 2n columns. Row operations are applied to the entire matrix. when you are ready to back substitute, divide the bottom row through by Mn,n. Subtract multiples of this row from the rows above to clear column n. Now the nth row and column of M look like the identity matrix. Move up to the previous row and divide through by Mn-1,n-1. Then subtract multiples of this row from previous rows to clear column n-1. Continue this process until M has become the identity matrix. At the same time, the identity matrix pasted onto the right of M has become R, the inverse of M.
Since M is nonsingular, each solution vector is unique, hence R is unique. Since MR = 1, write RMR = R. Now R has a unique right inverse; call it S. Multiply by S on the right to get RMRS = RS, or RM = 1. The right inverse is the left inverse. Nonsingular is synonymous with invertible.