code 33131;
procedure LINIGER2(X,XE,M,Y,SIGMA1,SIGMA2,F,EVALUATE,J,
               JACOBIAN,K,ITMAX,STEP,AETA,RETA,OUTPUT);
integer M,K,ITMAX;
real X,XE,SIGMA1,SIGMA2,STEP,AETA,RETA;
array Y,J;
boolean procedure EVALUATE;
real procedure F;
procedure JACOBIAN,OUTPUT;

begin integer I;
    real H,HL,B1,B2,P,Q,C0,C1,C2,C3,C4;
    boolean LAST;
    integer array PI[1:M];
    real array DY,YL,FL[1:M],A[1:M,1:M],AUX[1:3];

    procedure STEPSIZE;
    begin H:=STEP;
        if 1.1*H>=XE-X then 
        begin LAST:=true; H:=XE-X; X:=XE
        end else X:=X+H
    end STEPSIZE;

    procedure COEFFICIENT;
    begin real R1,R2,EX,ZETA,ETA,SINL,COSL,SINH,COSH,D;
        real procedure R(X); value X; real X;
        if X>40 then R:=X/(X-2) else 
        begin EX:=EXP(-X); R:=X*(1-EX)/(X-2+(X+2)*EX) end;

        B1:=H*SIGMA1;
        B2:=H*SIGMA2;
        if B1<.1 then begin P:=0; Q:=1/3; goto OUT end;
        if B2<0 then goto COMPLEX;
        if B1<1 or B2<.1 then goto THIRDORDER;
        if ABS(B1-B2)<B1*B1*"-6 then goto DOUBLEFIT;

        R1:=R(B1)*B1; R2:=R(B2)*B2;
        D:=B2*R1-B1*R2;
        P:=2*(R2-R1)/D;
        Q:=2*(B2-B1)/D;
THIRDORDER: Q:=1/3;
        P:=R(B1)/3-2/B1;
        goto OUT;
DOUBLEFIT: B1:=.5*(B1+B2);
        R1:=R(B1);
        if B1>40 then EX:=0;
        R2:=B1/(1-EX); R2:=1-EX*R2*R2;
        Q:=1/(R1*R1*R2);
        P:=R1*Q-2/B1;
        goto OUT;
COMPLEX: ETA:=ABS(B1*SIN(SIGMA2));
        ZETA:=ABS(B1*COS(SIGMA2));
        if ETA<B1*B1*"-6 then 
        begin B1:=B2:=ZETA; goto DOUBLEFIT end;
        if ZETA>40 then 
        begin P:=1-4*ZETA/B1/B1; Q:=4*(1-ZETA)/B1/B1+1 end else 
        begin EX:=EXP(ZETA);
            SINL:=SIN(ETA); COSL:=COS(ETA);
            SINH:=.5*(EX-1/EX); COSH:=.5*(EX+1/EX);
            D:=ETA*(COSH-COSL)-.5*B1*B1*SINL;
            P:=(ZETA*SINL+ETA*SINH-4*ZETA*ETA/B1/B1*(COSH-COSL))/D;
            Q:=ETA*((COSH-COSL-ZETA*SINH-ETA*SINL)*4/B1/B1+COSH+COSL)/D
        end;
OUT:    C0:=.25*H*H*(P+Q);
        C1:=.5*H*(1+P);
        C2:=H-C1;
        C3:=.25*H*H*(Q-P);
        C4:=.5*H*P;
        ELEMENTS OF MATRIX
    end COEFFICIENT;

    procedure ELEMENTS OF MATRIX;
    begin integer K;
        for I:=1 step 1 until M do 
        begin for K:=1 step 1 until M do 
            A[I,K]:=C0*MATMAT(1,M,I,K,J,J)-C1*J[I,K];
            A[I,I]:=A[I,I]+1
        end;
        AUX[2]:=0; DEC(A,M,AUX,PI)
    procedure NEWTON ITERATION;
    begin integer ITNUM; real JFL,ETA,DISCR;
        ITNUM:=0;
NEXT:   ITNUM:=ITNUM+1;
        if EVALUATE(ITNUM) then 
        begin JACOBIAN(J,Y); COEFFICIENT end 
        else if ITNUM=1 and H^=HL then COEFFICIENT;
        for I:=1 step 1 until M do FL[I]:=F(I);
        if ITNUM=1 then 
        begin for I:=1 step 1 until M do 
            begin JFL:=MATVEC(1,M,I,J,FL);
                DY[I]:=H*(FL[I]-C4*JFL);
                YL[I]:=Y[I]+C2*FL[I]+C3*JFL
            end 
        end else 
        for I:=1 step 1 until M do 
        DY[I]:=YL[I]-Y[I]+C1*FL[I]-C0*MATVEC(1,M,I,J,FL);
        SOL(A,M,PI,DY);
        for I:=1 step 1 until M do Y[I]:=Y[I]+DY[I];
        if ITNUM<ITMAX then 
        begin ETA:=SQRT(VECVEC(1,M,0,Y,Y))*RETA+AETA;
            DISCR:=SQRT(VECVEC(1,M,0,DY,DY));
            if ETA<DISCR then goto NEXT
        end 
    end NEWTON;

    LAST:=false; K:=0; HL:=0;
NEXT LEVEL:
    K:=K+1;
    STEPSIZE;
    NEWTON ITERATION;
    HL:=H;
    OUTPUT;
    if not LAST then goto NEXT LEVEL
end LINIGER2;
        eop