procedure GAUSSJ(VAR A:matrx2; N:integer; VAR B:array of real);
var
    IPIV,INDXR,INDXC:array[1..50] of integer;
    J,I,K,L,LL:integer;
    BIG,PIVINV,DUM:real; IROW,ICOL:integer;
begin
    For J:=1 To N do
      IPIV[J]:=0;
    For I:=1 To N do
    begin
        BIG:=0;
        For J:=1 To N do
        begin
            If IPIV[J] <> 1 Then
            begin
                For K:=1 To N do
                begin
                    If IPIV[K] = 0 Then
                    begin
                        If Abs(A[J, K]) >= BIG Then
                        begin
                            BIG:=Abs(A[J, K]);
                            IROW:=J;
                            ICOL:=K;
                        end;
                    end
                    Else if IPIV[K] > 1 Then
                        ShowMessage('Singular matrix.');
                end;
            end;
        end;
        IPIV[ICOL]:=IPIV[ICOL] + 1;
        If IROW <> ICOL Then
        begin
            For L:=1 To N do
            begin
                DUM:=A[IROW, L];
                A[IROW, L]:=A[ICOL, L];
                A[ICOL, L]:=DUM;
            end;
            DUM:=B[IROW];
            B[IROW]:=B[ICOL];
            B[ICOL]:=DUM;
        end;
        INDXR[I]:=IROW;
        INDXC[I]:=ICOL;
        If A[ICOL, ICOL] = 0 Then ShowMessage('Singular matrix.');
        PIVINV:=1 / A[ICOL, ICOL];
        A[ICOL, ICOL]:=1;
        For L:=1 To N do
          A[ICOL, L]:=A[ICOL, L] * PIVINV;
        B[ICOL]:=B[ICOL] * PIVINV;
        For LL:=1 To N do
        begin
            If LL <> ICOL Then
            begin
                DUM:=A[LL, ICOL];
                A[LL, ICOL]:=0;
                For L:=1 To N do
                    A[LL, L]:=A[LL, L] - A[ICOL, L] * DUM;
                B[LL]:=B[LL] - B[ICOL] * DUM;
            end;
        end;
    end;
    For L:=N DownTo 1 do
    begin
        If INDXR[L] <> INDXC[L] Then
        begin
            For K:=1 To N do
            begin
                DUM:=A[K, INDXR[L]];
                A[K, INDXR[L]]:=A[K, INDXC[L]];
                A[K, INDXC[L]]:=DUM;
            end;
        end;
    end;
end;