code 33120;
procedure EFERK(X,XE,M,Y,SIGMA,PHI,DERIVATIVE,J,JACOBIAN,
              K,L,AUT,AETA,RETA,HMIN,HMAX,LINEAR,OUTPUT);
value L; integer M,K,L;
real X,XE,SIGMA,PHI,AETA,RETA,HMIN,HMAX; array Y,J;
boolean AUT,LINEAR; procedure DERIVATIVE,JACOBIAN,OUTPUT;
begin integer M1,I;
    real H,B,B0,PHI0,COSPHI,SINPHI,ETA,DISCR,FAC,PI;
    boolean CHANGE,LAST;
    integer array P[1:L];
    real array BETA,BETHA[0:L],BETAC[0:L+3],K0,D,D1,D2[1:M],
        A[1:L,1:L],AUX[1:3];
    real procedure SUM(I,L,U,T); value L,U; integer I,L,U;
    real T;
    begin real S; S:=0;
        for I:=L step 1 until U do S:=S+T;
        SUM:=S
    end;
    procedure FORMBETA;
    if L=1 then 
    begin BETHA[1]:=(.5-(1-(1-EXP(-B))/B)/B)/B;
        BETA[1]:=(1/6-BETHA[1])/B
    end else 
    if L=2 then 
    begin real E,EMIN1; E:=EXP(-B); EMIN1:=E-1;
        BETHA[1]:=(1-(3+E+4*EMIN1/B)/B)/B;
        BETHA[2]:=(.5-(2+E+3*EMIN1/B)/B)/B/B;
        BETA[2]:=(1/6-BETHA[1])/B/B;
        BETA[1]:=(1/3-(1.5-(4+E+5*EMIN1/B)/B)/B)/B
    end else 
    begin real B0,B1,B2,A0,A1,A2,A3,C,D;
        BETAC[L-1]:=C:=D:=EXP(-B)/FAC;
        for I:=L-1 step -1 until 1 do 
        begin C:=I*B*C/(L-I); BETAC[I-1]:=D:=D*I+C end;
        B2:=.5-BETAC[2];
        B1:=(1-BETAC[1])*(L+1)/B;
        B0:=(1-BETAC[0])*(L+2)*(L+1)*.5/B/B;
        A3:=1/6-BETAC[3];
        A2:=B2*(L+1)/B;
        A1:=B1*(L+2)*.5/B;
        A0:=B0*(L+3)/3/B;
        D:=L/B;
        for I:=1 step 1 until L do 
        begin BETA[I]:=(A3/I-A2/(I+1)+A1/(I+2)-A0/(I+3))*D+BETAC[I+3];
            BETHA[I]:=(B2/I-B1/(I+1)+B0/(I+2))*D+BETAC[I+2];
            D:=D*(L-I)/I/B;
        end 
    procedure SOLUTIONOFCOMPLEXEQUATIONS;
    if L=2 then 
    begin real COS2PHI,COSA,SINA,E,ZI;
        PHI0:=PHI; COSPHI:=COS(PHI0); SINPHI:=SIN(PHI0);
        E:=EXP(B*COSPHI); ZI:=B*SINPHI-3*PHI0;
        SINA:=(if ABS(SINPHI)<"-6 then -E*(B+3)
              else E*SIN(ZI)/SINPHI);
        COS2PHI:=2*COSPHI*COSPHI-1;
        BETHA[2]:=(.5+(2*COSPHI+(1+2*COS2PHI+SINA)/B)/B)/B/B;
        SINA:=(if ABS(SINPHI)<"-6 then E*(B+4)
              else SINA*COSPHI-E*COS(ZI));
        BETHA[1]:=-(COSPHI+(1+2*COS2PHI+(4*COSPHI*COS2PHI+SINA)
                  /B)/B)/B;
        BETA[1]:=BETHA[2]+2*COSPHI*(BETHA[1]-1/6)/B;
        BETA[2]:=(1/6-BETHA[1])/B/B
    end else 

    begin integer J,C1;
        real C2,E,ZI,COSIPHI,SINIPHI,COSPHIL;
        real array D[1:L];
        procedure ELEMENTS OF MATRIX;
        begin PHI0:=PHI;
            COSPHI:=COS(PHI0); SINPHI:=SIN(PHI0);
            COSIPHI:=1; SINIPHI:=0;
            for I:=0 step 1 until L-1 do 
            begin C1:=4+I; C2:=1;
                for J:=L-1 step -2 until 1 do 
                begin  A[J,L-I]:=C2*COSIPHI;
                    A[J+1,L-I]:=C2*SINIPHI;
                    C2:=C2*C1; C1:=C1-1
                end;
                COSPHIL:=COSIPHI*COSPHI-SINIPHI*SINPHI;
                SINIPHI:=COSIPHI*SINPHI+SINIPHI*COSPHI;
                COSIPHI:=COSPHIL
            end;
            AUX[2]:=0; DEC(A,L,AUX,P)
        end EL OF MAT;
        procedure RIGHT HAND SIDE;
        begin E:=EXP(B*COSPHI);
            ZI:=B*SINPHI-4*PHI0;
            COSIPHI:=E*COS(ZI); SINIPHI:=E*SIN(ZI);
            ZI:=1/B/B/B;
            for J:=L step -2 until 2 do 
            begin D[J]:=ZI*SINIPHI;
                D[J-1]:=ZI*COSIPHI;
                COSPHIL:=COSIPHI*COSPHI-SINIPHI*SINPHI;
                SINIPHI:=COSIPHI*SINPHI+SINIPHI*COSPHI;
                COSIPHI:=COSPHIL; ZI:=ZI*B
            SINIPHI:=2*SINPHI*COSPHI;
            COSIPHI:=2*COSPHI*COSPHI-1;
            COSPHIL:=COSPHI*(2*COSIPHI-1);
            D[L]:=D[L]+SINPHI*(1/6+(COSPHI+(1+2*COSIPHI*(1+2*COSPHI/B))
                        /B)/B);
            D[L-1]:=D[L-1]-COSPHI/6-(.5*COSIPHI+(COSPHIL+(2*COSIPHI*
                        COSIPHI-1)/B)/B)/B;
            D[L-2]:=D[L-2]+SINPHI*(.5+(2*COSPHI+(2*COSIPHI+1)/B)/B);
            D[L-3]:=D[L-3]-.5*COSPHI-(COSIPHI+COSPHIL/B)/B;
            if L<5 then goto END;
            D[L-4]:=D[L-4]+SINPHI+SINIPHI/B;
            D[L-5]:=D[L-5]-COSPHI-COSIPHI/B;
            if L<7 then goto END;
            D[L-6]:=D[L-6]+SINPHI;
            D[L-7]:=D[L-7]-COSPHI;
        END:
        end RHS;
        if PHI0^=PHI then ELEMENTS OF MATRIX;
        RIGHT HAND SIDE;
        SOL(A,L,P,D);
        ZI:=1/B;
        for I:=1 step 1 until L do 
        begin BETA[I]:=D[L+1-I]*ZI;
            BETHA[I]:=(I+3)*BETA[I];
            ZI:=ZI/B
        end 
    end SOLOFEQCOM;

    procedure COEFFICIENT;
    begin B0:=B:=ABS(H*SIGMA);
        if B>=.1 then 
        begin if PHI^=PI and L=2 or ABS(PHI-PI)>.01 then 
            SOLUTION OF COMPLEX EQUATIONS else FORMBETA
        end else 
        begin for I:=1 step 1 until L do 
            begin BETHA[I]:=BETA[I-1];
                BETA[I]:=BETA[I-1]/(I+3);
            end 
        end 
    end COEFFICIENT;

    procedure LOCAL ERROR BOUND;
    ETA:=AETA+RETA*SQRT(VECVEC(1,M1,0,Y,Y));

    procedure STEPSIZE;
    begin LOCAL ERROR BOUND;
        if K=0 then 
        begin DISCR:=SQRT(VECVEC(1,M1,0,D,D)); H:=ETA/DISCR
        end else 
        begin DISCR:=H*SQRT(SUM(I,1,M1,(D[I]-D2[I])**2))/ETA;
            H:=H*(if LINEAR then 4/(4+DISCR)+.5
                              else 4/(3+DISCR)+1/3)
        if H<HMIN then H:=HMIN;
        if H>HMAX then H:=HMAX;
        B:=ABS(H*SIGMA);
        CHANGE:=ABS(1-B/B0)>.05 or PHI^=PHI0;
        if 1.1*H>=XE-X then 
        begin CHANGE:=LAST:=true; H:=XE-X end;
        if not CHANGE then H:=H*B0/B
    end STEPSIZE;

    procedure DIFFERENCE SCHEME;
    begin integer K;
        real BETAI,BETHAI;
        if M1<M then 
        begin D2[M]:=1; K0[M]:=Y[M]+2*H/3; Y[M]:=Y[M]+.25*H end;
        for K:=1 step 1 until M1 do 
        begin K0[K]:=Y[K]+2*H/3*D[K];
            Y[K]:=Y[K]+.25*H*D[K];
            D1[K]:=H*MATVEC(1,M,K,J,D);
            D2[K]:=D1[K]+D[K]
        end;
        for I:=0 step 1 until L do 
        begin BETAI:=4*BETA[I]/3; BETHAI:=BETHA[I];
            for K:=1 step 1 until M1 do D[K]:=H*D1[K];
            for K:=1 step 1 until M1 do 
            begin K0[K]:=K0[K]+BETAI*D[K];
                D1[K]:=MATVEC(1,M1,K,J,D);
                D2[K]:=D2[K]+BETHAI*D1[K]
            end 
        end;
        DERIVATIVE(K0);
        for K:=1 step 1 until M do Y[K]:=Y[K]+.75*H*K0[K]
    end DIFF SCHEME;

    B0:=PHI0:=-1; PI:=4*ARCTAN(1);
    BETAC[L]:=BETAC[L+1]:=BETAC[L+2]:=BETAC[L+3]:=0;
    BETA[0]:=1/6; BETHA[0]:=.5;
    FAC:=1; for I:=2 step 1 until L-1 do FAC:=I*FAC;
    M1:= if AUT then M else M-1;
    K:=0; LAST:=false;
 NEXT LEVEL:
    for I:=1 step 1 until M do D[I]:=Y[I];
    DERIVATIVE(D);
    if not LINEAR or K=0 then JACOBIAN(J,Y);
    STEPSIZE;
    if CHANGE then COEFFICIENT;
    OUTPUT;
    DIFFERENCE SCHEME;
    K:=K+1;
    X:=X+H;
    if not LAST then goto NEXT LEVEL;
END OF EFERK: OUTPUT;
end EFERK;
        eop