Questo è il mio primo programmino di una certa utilità (anche se non capita molto spesso di risolvere sistemi di equazioni), però io lo posto lo stesso: potrebbe sempre servire. E' stato realizzato in Java, ma si può facilmente adattare a qualsiasi linguaggio:
* @author Simone Salerno
 * @program Risolutore di sistemi di N equazioni lineari in N incognite con metodo di riduzione
//L'esempio si riferisce ad un sistema di 6 equazioni

public static void main(String[] args) {

        final int N = 6; //numero di equazioni (= incognite)
        int i, j;
        int m = -1, n = -1;
        final String[] Inc = {"X","Y","Z","T","K","H","M","N","O"};    //nomi incognite
        double[] inc = {0,0,0,0,0,0,0,0,0};                        //valori iniziali incognite
        final double[][] coeff = {{1,0,1,3,-1,-1,-5},{3,2,-1,-1,2,5,25},{2,1,2,3,10,-2,-10},
                              {10,-1,-1,-2,20,1,29},{20,10,3,-1,1,-1,9},{5,-5,-20,-20,1,3,40}}; 

//nell'ordine: {coeff. x, coeff.y, coeff.z, coeff.t, coeff. k, coeff. h, termine noto}
//SOLUZIONI: x=2;y=-3;z=1;t=-1;k=0;h=5
        
        while (m != N-1 && n != N-1) {
  
        //divide coeff. equazione per coeff. inc
        for (j=N-1; j>m; j--) {
            for (i=N; i>n; i--) 
                coeff [j][i] /= coeff [j][n+1];  }
        if (m - n == 0) m++;
        else n++;
        //sottrai coeff. "prima" equazione alle altre
        for (j=N-1; j>m; j--) {
            for (i=N; i>n; i--)
                coeff [j][i] -= coeff [j-1][i];        }
        if (m - n == 0) m++;
        else n++;
       }
        //Calcola soluzioni
        inc[N-1] = coeff [N-1][N] / coeff [N-1][N-1];
        for (i=N-2; i>-1; i--) {
            for (j=N, inc[i] = coeff[i][j]; j>0; j--) {
                if (i != j-1)
                inc[i] -= coeff [i][j-1]*inc[j-1];
            }
        }  
        //Assegna alle soluzioni non usate un valore di default
        for (i=N; i<inc.length; i++)
            inc[i] = -9999;
           
        //Stampa soluzioni
        for (i=0; i<Inc.length; i++)
        System.out.printf("%s = %3.1f; ", Inc[i], inc[i]);
    }
}
N. B.: i coefficienti delle x devono essere diversi da 0
N. B.: ci sarebbero altre restrizioni, ma per quanto è scemo il programma non perdo neanche tempo ad elencarle