code 33061;
procedure ARK (T, TE, M0, M, U, DERIVATIVE, DATA, OUT);
integer M0, M;
real T, TE;
array U, DATA;
procedure DERIVATIVE, OUT;

begin integer P, N, Q;
    own real EC0, EC1, EC2, TAU0, TAU1, TAU2, TAUS, T2;
    real THETANM1, TAU, BETAN, QINV, ETA;
    array MU, LAMBDA[1:DATA[1]], THETHA[0:DATA[1]], RO, R[M0:M];
    boolean START, STEP1, LAST;

    procedure INITIALIZE;
    begin integer I, J, K, L, N1; real S, THETA0;
        array ALFA[1:8, 1:DATA[1]+1], TH[1:8], AUX[1:3];

        real procedure LABDA(I, J); value I, J; integer I, J;
        LABDA:= if P < 3 then (if J =I-1 then MUI(I) else 0)
                else if P =3 then (if I =N then (if J=0
                then .25 else if J =N - 1 then .75
                else 0) else if J =0 then (if I =1
                then MUI(1) else .25) else if J =I - 1
                then LAMBDA[I] else 0) else 0;

        real procedure MUI(I); value I; integer I;
        MUI:= if I =N then 1 else 
              if I < 1 ! I > N then 0 else 
              if P < 3 then LAMBDA[I] else 
              if P =3 then LAMBDA[I] + .25 else 0;

        real procedure SUM(I, A, B, X);
        value B; integer I, A, B; real X;
        begin real S; S:= 0;
            for I:= A step 1 until B do S:= S + X;
            SUM:= S
        end SUM;
        N:= DATA[1]; P:= DATA[2]; EC1:= EC2 := 0;
        BETAN:= DATA[3];
        THETANM1:= if P=3 then .75 else 1;
        THETA0:= 1 - THETANM1; S:= 1;
        for J:= N - 1 step - 1 until 1 do 
        begin S:= - S * THETA0 + DATA[N + 10 - J];
            MU[J]:= DATA[N + 11 - J] / S;
            LAMBDA[J]:= MU[J] - THETA0
        end;
        for I:= 1 step 1 until 8 do 
        for J:= 0 step 1 until N do 
        ALFA[I, J + 1]:= if I = 1 then 1 else 
          if J = 0 then 0 else if I = 2 ! I = 4 ! I = 8 then 
          MUI(J) ** ENTIER((I + 2) / 3) else 
          if (I = 3 ! I = 6) & J > 1 then SUM(L, 1, J-1,
          LABDA(J, L) * MUI(L) ** ENTIER(I / 3)) else 
          if I = 5 & J > 2 then SUM(L, 2, J - 1, LABDA(J, L) *
          SUM(K, 1, L - 1, LABDA(L, K) * MUI(K))) else 
          if I = 7 & J > 1 then SUM(L, 1, J - 1, LABDA(J, L) *
          MUI(L)) * MUI(J) else 0;
        N1:=if N < 4 then N + 1 else if N < 7 then 4
          else 8;
        I:= 1;
        for S:= 1, .5, 1 / 6, 1 / 3, 1 / 24, 1 / 12, .125, .25 do 
        begin TH[I]:= S; I:= I + 1 end;
        if P = 3 & N < 7 then TH[1]:= TH[2]:= 0;
        AUX[2]:= " - 14; DECSOL(ALFA, N1, AUX, TH);
        INIVEC(0, N, THETHA, 0);
        DUPVEC(0, N1 - 1, 1, THETHA, TH);
        if ^ (P = 3 & N < 7) then 
        begin THETHA[0]:= THETHA[0] - THETA0;
            THETHA[N - 1]:= THETHA[N - 1] - THETANM1; Q:= P + 1
        end else Q:= 3;
        QINV:= 1 / Q;
        START:= DATA[8] = 0; DATA[10]:= 0; LAST:= false;
        DUPVEC(M0, M, 0, R, U); DERIVATIVE(T, R)
    end INITIALIZE

    procedure LOCAL ERROR CONSTRUCTION(I); value I; integer I;
    begin if THETHA[I] ^= 0 then 
        ELMVEC(M0, M, 0, RO, R, THETHA[I]);
        if I = N then 
        begin DATA[9]:= SQRT(VECVEC(M0, M, 0, RO, RO))* TAU;
            EC0:= EC1; EC1:= EC2; EC2:= DATA[9] / TAU ** Q
        end 
    end LEC;

    procedure STEPSIZE;
    begin real TAUACC, TAUSTAB, AA, BB, CC, EC;
        ETA:= SQRT(VECVEC(M0, M, 0, U, U)) * DATA[7] + DATA[6];
        if ETA > 0 then 
        begin if START then 
            begin if DATA[8] = 0 then 
                begin TAUACC:= DATA[5];
                    STEP1:= true 
                end else if STEP1 then 
                begin TAUACC:= (ETA / EC2) ** QINV;
                    if TAUACC > 10 * TAU2 then 
                    TAUACC:= 10 * TAU2 else STEP1:= false 
                end else 
                begin BB:= (EC2 - EC1) / TAU1; CC:= - BB * T2 + EC2;
                    EC:= BB * T + CC;
                    TAUACC:= if EC < 0 then TAU2 else 
                    (ETA / EC) ** QINV;
                    START:= false 
                end 
            end else 
            begin AA:= ((EC0 - EC1) / TAU0 + (EC2 - EC1) / TAU1)
                        / (TAU1 + TAU0);
                BB:= (EC2 - EC1) / TAU1 - (2 * T2 - TAU1) * AA;
                CC:= - (AA * T2 + BB) * T2 + EC2;
                EC:= (AA * T + BB) * T + CC;
                TAUACC:= if EC < 0 then 
                         TAUS else (ETA / EC) ** QINV;
                if TAUACC > 2 * TAUS then TAUACC:= 2 * TAUS;
                if TAUACC < TAUS / 2 then TAUACC:= TAUS / 2
            end 
        end else TAUACC:= DATA[5];
         if TAUACC < DATA[5] then TAUACC:= DATA[5];
        TAUSTAB:= BETAN / DATA[4]; if TAUSTAB < DATA[5] then 
        begin DATA[10]:= 1; goto ENDARK end;
        TAU:= if TAUACC > TAUSTAB then TAUSTAB else TAUACC;
        TAUS:= TAU; if TAU >= TE - T then 
        begin TAU:= TE - T; LAST:= true end;
        TAU0:= TAU1; TAU1:= TAU2; TAU2:= TAU
    end STEPSIZE

    procedure DIFFERENCE SCHEME;
    begin integer I, J;
        real MT, LT;
        MULVEC(M0, M, 0, RO, R, THETHA[0]);
        if P = 3 then ELMVEC(M0, M, 0, U, R, .25 * TAU);
        for I:= 1 step 1 until N - 1 do 
        begin MT:= MU[I] * TAU; LT:= LAMBDA[I] * TAU;
            for J:= M0 step 1 until M do 
            R[J]:= LT * R[J] + U[J];
            DERIVATIVE(T + MT, R); LOCAL ERROR CONSTRUCTION(I)
        end;
        ELMVEC(M0, M, 0, U, R, THETANM1 * TAU);
        DUPVEC(M0, M, 0, R, U); DERIVATIVE(T + TAU, R);
        LOCAL ERROR CONSTRUCTION(N); T2:= T;
        if LAST then 
        begin LAST:= false; T:= TE end else T:= T + TAU;
        DATA[8]:= DATA[8]+1
    end DIFSCH;

    INITIALIZE;

  NEXT STEP:
    STEPSIZE; DIFFERENCE SCHEME; OUT;
    if T ^= TE then goto NEXT STEP;

  ENDARK:
end ARK;
        eop