procedure ADI(A, B, C, D, E, F, G:matrx2;var U:matrx2;
                           JMAX, K:integer; ALPHA, BETA, EPS:real);
const
    JJ = 100;    KK = 6;   MAXITS = 100;  ZERO = 0;
    TWO = 2;     HALF = 0.5;
var
    AA, BB, CC, RR, UU:array[0..100] of real;  PSI,S:matrx2;
    ALPH, BET:array[0..6] of real; R:array[0..32] of real;
    J,L,N,NR,NRR,K1,NITS,KITS,NEXT1:integer;
    AB,ANORMG,DISC,RFACT,ANORM,RESID,ANROM:real;
begin
    SetLength(PSI,101,101);
    SetLength(S,33,7);
    NRR:= Trunc(EXP((KK - 1)*Ln(2)));
    If JMAX > JJ Then ShowMessage('Increase JJ');
    If K > KK - 1 Then ShowMessage('Increase KK');
    K1:=K + 1;
    NR:= Trunc(Exp(K*Ln(2)));
    ALPH[1]:=ALPHA;
    BET[1]:=BETA;
    For J:=1 To K do
    begin
        ALPH[J + 1]:=Sqrt(ALPH[J] * BET[J]);
        BET[J + 1]:=HALF * (ALPH[J] + BET[J]);
    end;
    S[1, 1]:=Sqrt(ALPH[K1] * BET[K1]);
    For J:=1 To K do
    begin
        AB:=ALPH[K1 - J] * BET[K1 - J];
        For N:=1 To Trunc(Exp((J - 1)*Ln(2))) do
        begin
            DISC:=Sqrt(Sqr(S[N, J]) - AB);
            S[2 * N, J + 1]:=S[N, J] + DISC;
            S[2 * N - 1, J + 1]:=AB / S[2 * N, J + 1];
        end;
    end;
    For N:=1 To NR do
        R[N]:=S[N, K1];
    ANORMG:=ZERO;
    For J:=2 To JMAX - 1 do
    begin
      For L:=2 To JMAX - 1 do
      begin
        ANORMG:=ANORMG + AbS(G[J, L]);
        PSI[J, L]:=-D[J, L] * U[J, L - 1] + (R[1] - E[J, L]) * U[J, L];
        PSI[J, L]:=PSI[J, L] - F[J, L] * U[J, L + 1];
      end;
    end;
    NITS:=MAXITS div NR;
    For KITS:=1 To NITS do
    begin
        For N:=1 To NR do
        begin
            If N = NR Then
                NEXT1:=1
            Else
                NEXT1:=N + 1;
            RFACT:=R[N] + R[NEXT1];
            For L:=2 To JMAX - 1 do
            begin
                For J:=2 To JMAX - 1 do
                begin
                    AA[J - 1]:=A[J, L];
                    BB[J - 1]:=B[J, L] + R[N];
                    CC[J - 1]:=C[J, L];
                    RR[J - 1]:=PSI[J, L] - G[J, L];
                end;
                TRIDAG(AA, BB, CC, RR, UU, JMAX - 2);
                For J:=2 To JMAX - 1 do
                    PSI[J, L]:=-PSI[J, L] + TWO * R[N] * UU[J - 1];
            end;
            For J:=2 To JMAX - 1 do
            begin
                For L:=2 To JMAX - 1 do
                begin
                    AA[L - 1]:=D[J, L];
                    BB[L - 1]:=E[J, L] + R[N];
                    CC[L - 1]:=F[J, L];
                    RR[L - 1]:=PSI[J, L];
                end;
                TRIDAG(AA, BB, CC, RR, UU, JMAX - 2);
                For L:=2 To JMAX - 1 do
                begin
                    U[J, L]:=UU[L - 1];
                    PSI[J, L]:=-PSI[J, L] + RFACT * UU[L - 1];
                end;
            end;
        end;
        ANORM:=ZERO;
        For J:=2 To JMAX - 1 do
        begin
          For L:=2 To JMAX - 1 do
          begin
            RESID:=A[J, L] * U[J - 1, L] + (B[J, L] + E[J, L]) * U[J, L];
            RESID:=RESID + C[J, L] * U[J + 1, L];
            RESID:=RESID + D[J, L] * U[J, L - 1] * U[J, L - 1];
            RESID:=RESID + F[J, L] * U[J, L + 1] + G[J, L];
            ANROM:=ANORM + AbS(RESID);
          end;
        end;
        If ANORM < EPS * ANORMG Then Exit;
    end;
    ShowMessage(' MAXITS exceeded');
end;