code 33160;
    procedure EFSIRK(X, XE, M, Y, DELTA, DERIVATIVE, JACOBIAN, J,
                       N, AETA, RETA, HMIN, HMAX, LINEAR, OUTPUT);
    value M; integer M, N;
    real X, XE, DELTA, AETA, RETA, HMIN, HMAX;
    procedure DERIVATIVE, JACOBIAN, OUTPUT;
    boolean LINEAR;
    array Y, J;

    begin integer K, L;
        real STEP, H, MU0, MU1, MU2, THETA0, THETA1, NU1, NU2,
        NU3, YK, FK, C1, C2, D;
        array F, K0, LABDA[1 : M], J1[1 : M,1 : M], AUX[1 : 7];
        integer array RI, CI[1 : M];
        boolean LIN;
        real procedure STEPSIZE;
        begin real DISCR, ETA, S;
            if LINEAR then S:= H:= HMAX else 
            if N = 1 or HMIN = HMAX then S:= H:= HMIN else 
            begin ETA:= AETA + RETA * SQRT(VECVEC(1, M, 0, Y, Y));
                C1:= NU3 * STEP; for K:= 1 step 1 until M do 
                LABDA[K]:= LABDA[K] + C1 * F[K] - Y[K];
                DISCR:= SQRT(VECVEC(1, M, 0, LABDA, LABDA));
                S:= H:= (ETA / (0.75 * (ETA + DISCR)) + 0.33) * H;
                if H < HMIN then S:= H:= HMIN else 
                if H > HMAX then S:= H:= HMAX
            end;
            if X + S > XE then S:= XE - X;
            LIN:= STEP = S and LINEAR; STEPSIZE:= S
        end STEPSIZE;

        procedure COEFFICIENT;
        begin real Z1, E, ALPHA1, A, B;
            own real Z2;
            Z1:= STEP * DELTA; if N = 1 then Z2:= Z1 + Z1;
            if ABS(Z2 - Z1) > " - 6 * ABS(Z1) or Z2 > - 1 then 
            begin A:= Z1 * Z1 + 12; B:= 6 * Z1;
                if ABS(Z1) < 0.1 then 
                ALPHA1:= (Z1 * Z1 / 140 - 1) * Z1 / 30 else 
                if Z1 < - "14 then ALPHA1:= 1 / 3 else 
                if Z1 < - 33 then 
                ALPHA1:= (A + B) / (3 * Z1 * (2 + Z1)) else 
                begin E:= if Z1 < 230 then EXP(Z1) else "100;
                    ALPHA1:= ((A - B) * E - A - B) /
                             (((2 - Z1) * E - 2 - Z1) * 3 * Z1)
                MU2:= (1 / 3 + ALPHA1) * 0.25;
                MU1:= - (1 + ALPHA1) * 0.5;
                MU0:= (6 * MU1 + 2) / 9; THETA0:= 0.25;
                THETA1:= 0.75; A:= 3 * ALPHA1;
                NU3:= (1 + A) / (5 - A) * 0.5; A:= NU3 + NU3;
                NU1:= 0.5 - A; NU2:= (1 + A) * 0.75;
                Z2:= Z1
            end 
        end COEFFICIENT;

        procedure DIFFERENCE SCHEME;
        begin DERIVATIVE(F); STEP:= STEPSIZE;
            if not LINEAR or N = 1 then JACOBIAN(J, Y);
            if not LIN then 
            begin COEFFICIENT;
                C1:= STEP * MU1; D:= STEP * STEP * MU2;
                for K:= 1 step 1 until M do 
                begin for L:= 1 step 1 until M do 
                    J1[K,L]:= D * MATMAT(1, M, K, L, J, J) +
                              C1 * J[K,L];
                    J1[K,K]:= J1[K,K] + 1
                end;
                GSSELM(J1, M, AUX, RI, CI)
            end;
            C1:= STEP * STEP * MU0; D:= STEP * 2 / 3;
            for K:= 1 step 1 until M do 
            begin K0[K]:= FK:= F[K];
                LABDA[K]:= D * FK + C1 * MATVEC(1, M, K, J, F)
            end;
            SOLELM(J1, M, RI, CI, LABDA);
            for K:= 1 step 1 until M do F[K]:= Y[K] + LABDA[K];
            DERIVATIVE(F);
            C1:= THETA0 * STEP; C2:= THETA1 * STEP; D:= NU1 * STEP;
            for K:= 1 step 1 until M do 
            begin YK:= Y[K]; FK:= F[K];
                LABDA[K]:= YK + D * FK + NU2 * LABDA[K];
                Y[K]:= F[K]:= YK + C1 * K0[K] + C2 * FK
            end 
        end DIFFERENCE SCHEME;

        AUX[2]:= "-14; AUX[4]:= 8;
        for K:= 1 step 1 until M do F[K]:= Y[K];
        N:= 0; OUTPUT; STEP:= 0;
    NEXT STEP: N:= N + 1;
        DIFFERENCE SCHEME; X:= X + STEP; OUTPUT;
        if X < XE then goto NEXT STEP
    end EFSIRK;
        eop