code 33070;
procedure EFRK(T,TE,M0,M,U,SIGMA,PHI,DIAMETER,DERIVATIVE,K,STEP,R,L,
        BETA,THIRDORDER,TOL,OUTPUT);
value R,L;
integer M0,M,K,R,L;
real T,TE,SIGMA,PHI,DIAMETER,STEP,TOL;
array U,BETA;
boolean THIRDORDER;
procedure DERIVATIVE,OUTPUT;
begin integer N;
    real THETA0,THETANM1,H,B,B0,PHI0,PHIL,PI,COSPHI,SINPHI,EPS,BETAR;
    boolean FIRST,LAST,COMPLEX,CHANGE;
    integer array P[1:L];
    real array MU,LABDA[0:R+L-1],PT[0:R],FAC,BETAC[0:L-1],RL[M0:M],
                   A[1:L,1:L],AUX[0:3];
    procedure FORM CONSTANTS;
    begin integer I;
        FIRST:=false;
        FAC[0]:=1;
        for I:=1 step 1 until L-1 do FAC[I]:=I*FAC[I-1];
        PT[R]:=L*FAC[L-1];
        for I:=1 step 1 until R do 
        PT[R-I]:=PT[R-I+1]*(L+I)/I
    end FORM CONSTANTS;

    procedure FORM BETA;
    begin integer I,J; real BB,C,D;
        if FIRST then FORM CONSTANTS;
        if L=1 then 
        begin C:=1-EXP(-B);
            for J:=1 step 1 until R do C:=BETA[J]-C/B;
            BETA[R+1]:=C/B
        end else 
        if B>40 then 
        begin for I:=R+1 step 1 until R+L do 
            begin C:=0;
                for J:=0 step 1 until R do 
                C:=BETA[J]*PT[J]/(I-J)-C/B;
                BETA[I]:=C/B/FAC[L+R-I]/FAC[I-R-1]
            end;
        end else 
        begin D:=C:=EXP(-B); BETAC[L-1]:=D/FAC[L-1];
            for I:=1 step 1 until L-1 do 
            begin C:=B*C/I; D:=D+C; BETAC[L-1-I]:=D/FAC[L-1-I] end;
            BB:=1;
            for I:=R+1 step 1 until R+L do 
            begin C:=0;
                for J:=0 step 1 until R do 
                C:=(BETA[J]-(if J<L then BETAC[J] else 0))*
                   PT[J]/(I-J)-C/B;
                BETA[I]:=C/B/FAC[L+R-I]/FAC[I-R-1]+
                         (if I<L then BB*BETAC[I] else 0);
                BB:=BB*B
            end 
        end 
    end FORM BETA;

    procedure SOLUTION OF COMPLEX EQUATIONS;
    begin integer I,J,C1,C3;
        real C2,E,B1,ZI,COSIPHI,SINIPHI,COSPHIL;
        real array D[1:L];

        procedure ELEMENTS OF MATRIX;
        begin PHIL:=PHI0;
            COSPHI:=COS(PHIL); SINPHI:=SIN(PHIL);
            COSIPHI:=1; SINIPHI:=0;
            for I:=0 step 1 until L-1 do 
            begin C1:=R+1+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:=C1*C2; 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 RIGHTHANDSIDE;
        begin E:=EXP(B*COSPHI);
            B1:=B*SINPHI-(R+1)*PHIL;
            COSIPHI:=E*COS(B1); SINIPHI:=E*SIN(B1);
            B1:=1/B; ZI:=B1**R;
            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
            end;
            COSIPHI:=ZI:=1; SINIPHI:=0;
            for I:=R step -1 until 0 do 
            begin C1:=I; C2:=BETA[I];
                C3:=if 2*I>L-2 then 2 else L-2*I;
                COSPHIL :=COSIPHI*COSPHI-SINIPHI*SINPHI;
                SINIPHI:=COSIPHI*SINPHI+SINIPHI*COSPHI;
                COSIPHI:=COSPHIL;
                for J:=L step -2 until C3 do 
                begin D[J]:=D[J]+ZI*C2*SINIPHI;
                    D[J-1]:=D[J-1]-ZI*C2*COSIPHI;
                    C2:=C2*C1; C1:=C1-1
                end;
                ZI:=ZI*B1
            end 
        end RIGHT HAND SIDE;

        if PHI0^=PHIL then ELEMENTS OF MATRIX;
        RIGHTHANDSIDE;
        SOL(A,L,P,D);
        for I:=1 step 1 until L do BETA[R+I]:=D[L+1-I]*B1
    procedure COEFFICIENT;
    begin integer J,K; real C;
        B0:=B; PHI0:=PHI;
        if B>=1 then 
        begin if COMPLEX then SOLUTION OF COMPLEX EQUATIONS
            else FORM BETA
        end;
        LABDA[0]:=MU[0]:=0;
        if THIRDORDER then 
        begin THETA0:=.25; THETANM1:=.75;
            if B<1 then 
            begin C:=MU[N-1]:=2/3; LABDA[N-1]:=5/12;
                for J:=N-2 step -1 until 1 do 
                begin C:=MU[J]:=C/(C-.25)/(N-J+1);
                    LABDA[J]:=C-.25
                end 
            end else 
            begin C:=MU[N-1]:=BETA[2]*4/3; LABDA[N-1]:=C-.25;
                for J:=N-2 step -1 until 1 do 
                begin C:=MU[J]:=C/(C-.25)*BETA[N-J+1]/BETA[N-J]/
                           (if J<L then B else 1);
                    LABDA[J]:=C-.25
                end 
            end 
        end else 
        begin THETA0:=0; THETANM1:=1;
            if B<1 then 
            begin for J:=N-1 step -1 until 1 do 
                MU[J]:=LABDA[J]:=1/(N-J+1)
            end else 
            begin LABDA[N-1]:=MU[N-1]:=BETA[2];
                for J:=N-2 step -1 until 1 do 
                MU[J]:=LABDA[J]:=BETA[N-J+1]/BETA[N-J]/
                       (if J<L then B else 1)
            end 
        end 
    procedure STEPSIZE;
    begin real D,HSTAB,HSTABINT;
        H:=STEP;
        D:=ABS(SIGMA*SIN(PHI));
        COMPLEX:=L//2*2=L and 2*D>DIAMETER;
        if DIAMETER>0 then 
        HSTAB:=(SIGMA**2/(DIAMETER*(DIAMETER*.25+D)))**(L*.5/R)/
               BETAR/SIGMA
        else HSTAB:=H;
        D:= if THIRDORDER then (2*TOL/EPS/BETA[R])**(1/(N-1))*
            4**((L-1)/(N-1)) else (TOL/EPS)**(1/R)/BETAR;
        HSTABINT:= ABS(D/SIGMA);
        if H>HSTAB then H:=HSTAB;
        if H>HSTABINT then H:=HSTABINT;
        if T+H>TE*(1-K*EPS) then 
        begin LAST:=true; H:=TE-T end;
        B:=H*SIGMA; D:=DIAMETER*.1*H; D:=D*D;
        if H<T*EPS  then goto ENDOFEFRK;
        CHANGE:=B0=-1 or ((B-B0)*(B-B0)+B*B0*(PHI-PHI0)*(PHI-PHI0)>D)
    end STEPSIZE;

    procedure DIFFERENCESCHEME ;
    begin integer I,J; real MT,LT,THT;
        I:=-1;
        NEXTTERM:
        I:=I+1; MT:=MU[I]*H; LT:=LABDA[I]*H;
        for J:=M0 step 1 until M do RL[J]:=U[J]+LT*RL[J];
        DERIVATIVE(T+MT,RL);
        if I=0 or I=N-1 then 
        begin THT:=if I=0 then THETA0*H else THETANM1*H;
            ELMVEC(M0,M,0,U,RL,THT)
        end;
        if I<N-1 then goto NEXTTERM;
        T:=T+H
    end DIFFERENCE SCHEME;

    N:=R+L; FIRST:=true; B0:=-1; BETAR:=BETA[R]**(1/R);
    LAST:=false; EPS:=2.0**(-48); PI:=PHI0:=PHIL:=4*ARCTAN(1);
    INIVEC(M0, M, RL, 0);
NEXTLEVEL:
    STEPSIZE;
    if CHANGE then COEFFICIENT;
    K:=K+1;
    DIFFERENCE SCHEME;
    OUTPUT;
    if not LAST then goto NEXTLEVEL;
ENDOFEFRK:
end EXPONENTIALLY FITTED RUNGE KUTTA;
        eop