Here is presented the so-called Gauss's method to solve the system of linear equations. You can download program. The algorithm have been realized on C.

Assume we have the system of N linear equations

a₁₁x₁ + a₁₂x₂ + a₁₃x₃ + ... a_1Nx_N = b₁
a₂₁x₁ + a₂₂x₂ + a₂₃x₃ + ... a_2Nx_N = b₂
a₃₁x₁ + a₃₂x₂ + a₃₃x₃ + ... a_3Nx_N = b₃
...
a_N1x₁ + a_N2x₂ + a_N3x₃ + ... a_NNx_N = b_N

where x_i - variable, a_ij - are coefficients before variable, b_i - are free terms, i,j are running from 1 to N.

The aim is to find x_i on basis of a_ij and b_i.

The Gauss's method consists in the reduction of system to simplest triangle form. The triangle form of system looks like

a11x1 +	a12x2 +	a13x3 +	...	a1NxN = b1
	a22x2 +	a23x3 +	...	a2NxN = b2
		a33x3 +	...	a3NxN = b3
...
			...	aNNxN = bN

The main feature of this system is that all a_ij are zeroes where j<i

The solution of triangle form of system is simple. We can easily find x_N= b_N / a_NN from the last equation. After substitution of x_N into the next to the last equation we can find x_N-1. And so on, and so on until we find x₁. This is the method of backward-sweep.

Now consider the method to transform the system into triangle form.

As the linear algebra teachs (English reference required) it's known that we can add to one line of system the any linear combination of any other lines of system. And the solution of system will not changed. The linear combination is the sum of lines multipied to some value.

If we want the second equation where the first term near x₁ is zero we can add to second line the first line multiplyed to some value M.

(a₁₁x₁ + a₁₂x₂ + a₁₃x₃ + ... a_1Nx_N = b₁)*M +
a₂₁x₁ + a₂₂x₂ + a₂₃x₃ + ... a_2Nx_N = b₂

Here we have

(a₁₁*л + a₂₁) x₁ + ... = b₁*M + b₂

To make term near x₁ to be equal zero we have to assign M = - a₂₁ / a₁₁. After the addition we substitute the result equation instead of second one. The third equation can be similarly transformed after the addition of the first equation multipyed by M = - a₃₁ / a₁₁. After that we have zero instead of term near x₁. After the operation described above have been applyed to all other equations we will have the next system:

a11x1 +	a12x2 +	a13x3 +	...	a1NxN = b1
	a22x2 +	a23x3 +	...	a2NxN = b2
	a32x2 +	a33x3 +	...	a3NxN = b3
...
	aN2x2 +	aN3x3 +	...	aNNxN = bN

After that we will make terms near x₂ to be zero in 3rd, 4th, N-th equation. To make it we have to add 2nd equation multiplyed on M = - a_j2 / a₂₂ to j-th equation. After that stage we will have the system without terms near x₂ at lines with numbers more than 2.

And so on... We will make this procedure for 3rd, 4th term untill the system will end.

The text of program.

//HOWTO solve system of linear equations
//ax=b
//where a[N][N] is square matrix,
//b[N] and x[N] are columns

//data are stored in input text file:

//input file format:
//a_11 a_12 a_13... a_1N b_1
//a_21 a_22 a_23... a_2N b_2
//...
//a_N1 a_N2 a_N3... a_NN b_N

#include <stdio.h>
#include <process.h>
float **a, *b, *x;
int N;
char filename[256];
FILE* InFile=NULL;
void count_num_lines(){
   //count number of lines in input file - number of equations
   int nelf=0;       //non empty line flag
   do{
       nelf = 0;
       while(fgetc(InFile)!='\n' && !feof(InFile)) nelf=1;
       if(nelf) N++;
   }while(!feof(InFile));
}

void freematrix(){
   //free memory for matrixes
   int i;
   for(i=0; i<N; i++){
       delete [] a[i];
   }
   delete [] a;
   delete [] b;
   delete [] x;
}

void allocmatrix(){
   //allocate memory for matrixes
   int i,j;
   x = new float[N];
   b = new float[N];
   a = new float*[N];
   if(x==NULL || b==NULL || a==NULL){
       printf("\nNot enough memory to allocate for %d equations.\n", N);
       exit(-1);
   }
   for(i=0; i<N; i++){
       a[i] = new float[N];
       if(a[i]==NULL){
	   printf("\nNot enough memory to allocate for %d equations.\n", N);
       }
   }
   for(i=0; i<N; i++){
       for(j=0; j<N; j++){
	   a[i][j]=0;
       }
       b[i]=0;
       x[i]=0;
   }
}

void readmatrix(){
   int i=0,j=0;
   //read matrixes a and b from input file
   for(i=0; i<N; i++){
       for(j=0; j<N; j++){
	   fscanf(InFile, "%f", &a[i][j]);
       }
       fscanf(InFile, "%f", &b[i]);
   }
}

void printmatrix(){
   //print matrix "a"
   int i=0,j=0;
   printf("\n");
   for(i=0; i<N; i++){
       for(j=0; j<N; j++){
	   printf(" %+f*X%d", a[i][j], j);
       }
       printf(" =%f\n", b[i]);
   }
}

void testsolve(){
   //test that ax=b
   int i=0,j=0;
   printf("\n");
   for(i=0; i<N; i++){
       float s = 0;
       for(j=0; j<N; j++){
	   s += a[i][j]*x[j];
       }
       printf("%f\t%f\n", s, b[i]);
   }
}
void printresult(){
   int i=0;
   printf("\n");
   printf("Result\n");
   for(i=0; i<N; i++){
       printf("X%d = %f\n", i, x[i]);
   }
}
void diagonal(){
   int i, j, k;
   float temp=0;
   for(i=0; i<N; i++){
       if(a[i][i]==0){
	   for(j=0; j<N; j++){
	       if(j==i) continue;
	       if(a[j][i] !=0 && a[i][j]!=0){
		   for(k=0; k<N; k++){
		       temp = a[j][k];
		       a[j][k] = a[i][k];
		       a[i][k] = temp;
		   }
		   temp = b[j];
		   b[j] = b[i];
		   b[i] = temp;
		   break;
	       }
	   }
       }
   }
}
void cls(){
   for(int i=0; i<25; i++) printf("\n");
}
void main(){
   int i=0,j=0, k=0;
   cls();
   do{
       printf("\nInput filename: ");
       scanf("%s", filename);
       InFile = fopen(filename, "rt");
   }while(InFile==NULL);
   count_num_lines();
   allocmatrix();
   rewind(InFile);
   //read data from file
   readmatrix();
   //check if there are 0 on main diagonal and exchange rows in that case
   diagonal();
   fclose(InFile);
   printmatrix();
   //process rows
   for(k=0; k<N; k++){
       for(i=k+1; i<N; i++){
	   if(a[k][k]==0){
	       printf("\nSolution is not exist.\n");
	       return;
	   }
	   float M = a[i][k] / a[k][k];
	   for(j=k; j<N; j++){
	       a[i][j] -= M * a[k][j];
	   }
	   b[i] -= M*b[k];
       }
   }
   printmatrix();
   for(i=N-1; i>=0; i--){
       float s = 0;
       for(j = i; j<N; j++){
	   s = s+a[i][j]*x[j];
       }
       x[i] = (b[i] - s) / a[i][i];
   }

   InFile = fopen(filename, "rt");
   readmatrix();
   fclose(InFile);
   printmatrix();
   testsolve();
   printresult();
   freematrix();
}

Thus the method have been described above I will describe the interaction of user and program.

Run the compiled executive file to start working with program. The program will ask you where to take the coefficients of equations. You can write text file (use notepad, not Word) with the coefficients separated by spaces in lines:

a11	a12	a13	...	a1N	b1
a21	a22	a23	...	a2N	b2
a31	a32	a33	...	a3N	b3
aN1	aN2	aN3	...	aNN	bN

For example:

2	-3	1	5
-1	-1	2	0
1	2	-1	3

You have to save file in the same directory where the progam exists. In other case the program will not find it. After the working the program will print result like:

X0 = 3.000000
X1 = 1.000000
X2 = 2.000000

It is the solution of system e.g. the set of unknown у.

I'll shortly describe the purpose of subroutines of program:

• void count_num_lines() - to calculate the number of equations in system
• void allocmatrix() - to get the memory for arrays for store the coefficients of equations, the free terms and the solution.
• void readmatrix() - to read the coefficients and the free terms from file
• void printmatrix() - to print the system
• void diagonal() - to avoid the zeroes on the main diagonal to escape the possibility of division on zero
• void testsolve() - to substitute the solution into system and to print the resulting right part of equation and the initial free terms. The two columns have to be equal in the case of good solution.
• void printresult() - to print the result column of solution
• void freematrix() - to free the memory allocated before
• void cls() - to clear the screen in the begining of program work
• void main() - the main function which call all of the subroutines described above, which transforms the system to triangle form and makes the backward-sweep.

The program have comments which ease it's understanding.

I hope you have understand everything.