Resolver sistemas lineares

Program Sistema_Linear;

    type vetor = array[1..20] of extended;
        matriz = array[1..20,1..20] of extended;

    var Constantes: vetor;
        Mampliada: matriz;
        Mauxiliar: matriz;
        Cons_auxiliar: vetor;
        i, j, N, Num_escolha: integer;
        Escolha: char;

function compara_vetor (Vetora, Vetorb: vetor): extended;
    var Subtrai: vetor;
        Max: extended;
    begin
        max := 0;
        for i := 1 to N do
            begin
                Subtrai[i] := Vetora[i] - Vetorb[i];
            end;
        for i := 1 to N do
            begin
                if abs(Subtrai[i]) > Max then
                    begin
                        Max := abs(Subtrai[i]);
                    end;
                compara_vetor := max;
            end;
    end;

procedure escalona_sup (MatrizA: matriz; Vetorconst: vetor);

    var i, j, k: integer;
        aux: extended;
    begin
        for k := 1 to (N - 1) do
            begin
                for j := (k + 1) to N do
                    begin
                        aux := MatrizA[k, j]/MatrizA[j, j];
                        for j := (k + 1) to N do
                            begin
                                MatrizA[k, i] := MatrizA[k, i] - aux * MatrizA[j, i];
                            end;
                        Vetorconst[k] := Vetorconst[k] - aux * Vetorconst[j];
                    end;
            end;
    end;

procedure escalona_inf (MatrizA: matriz; Vetorconst: vetor);
    var i, k : integer;
        s: extended;
    begin
        for k := 1 to (N - 1) do
            begin
                for i := (k + 1) to N do
                    begin
                        s := MatrizA[i, k]/MatrizA[k, k];
                        MatrizA[i, k] := 0;
                        for j := (k + 1) to N do
                            begin
                                MatrizA[i, j] := MatrizA[i, j] - s * MatrizA[k, j];
                            end;
                        Vetorconst[i] := Vetorconst[i] - s * Vetorconst[k];
                    end;
            end;
    end;

procedure gauss (MatrizAmp: matriz; Constante: vetor);
    var i, k: integer;
        Resposta: vetor;
        aux: extended;
    begin
        writeln(' ');
        writeln('Gauss');
        Mauxiliar := MatrizAmp;
        Cons_auxiliar := Constante;
        escalona_inf (Mauxiliar, Cons_auxiliar);
        Resposta[N] := Cons_auxiliar[N]/Mauxiliar[N, N];
        for i := (N - 1) downto 1 do
            begin
                aux := 0;
                for j := (i + 1) to N do
                    begin
                        aux := aux + Mauxiliar[i, j] * Resposta[j];
                        Resposta[i] := (Mauxiliar[i] - aux)/Mauxiliar[i, i];
                    end;
                for i := 1 to N do
                    begin
                        writeln('X', i, ' = ', Resposta[i]);
                        writeln(' ');
                    end;
            end;
    end;

begin
    repeat
    clrscr;
    write('Quantidade de variáveis = ');
    readln(N);
    for i := 1 to N do
        begin
            writeln(' ');
            writeln('Fornecer os coeficientes da equacao ', i, ':');
            for j := 1 to N do
                begin
                    write('Coeficiente de X', j, ' = ');
                    readln(Mampliada[i,j]);
                end;
            write('Valor da constante da equacao ', i,' = ');
            readln(Constantes[i]);
    end;
    writeln(' ');
    writeln('Escolha o metodo de resolucao:');
    writeln(' ');
    writeln('1 - Eliminacao de Gauss');
    writeln('2 - Metodo de eliminacao de Gauss-Jordan');
    writeln('3 - Metodo de Jacobi');
    writeln('4 - Metodo de Gauss-Seidel');
    writeln(' ');
    repeat
    write(' Digite o numero do metodo escolhido (1-4): ');
    readln(Num_escolha);
    case num of
        1: gauss(Mampliada, Constantes);
        2: jordan(Mampliada, Constantes);
        3: jacobi(Mampliada, Constantes);
        4: seidel(Mampliada, Constantes);
    end;
    write('Escolher outro metodo? (S/N): ');
    readln(Escolha);
    Escolha := upcase(Escolha);
    until Escolha = 'N';
    write('Sair do programa? (S/N): ');
    readln(Escolha);
    resposta := upcase(Escolha);
    until Escolha = 'S';
end.