code 33135;
procedure IMPEX (N, T0, TEND, Y0, DERIV, AVAILABLE, H0, HMAX,
                  PRESCH, EPS, WEIGHTS, UPDATE, FAIL, CONTROL);
value N;
integer N;
real T0,TEND,H0,HMAX,EPS;
boolean PRESCH,FAIL;
array Y0,WEIGHTS;
boolean procedure AVAILABLE;
procedure DERIV,UPDATE,CONTROL;
begin integer I,K,ECI;
    real T,T1,T2,T3,TP,H,H2,HNEW,ALF,LQ;
    array Y,Z,S1,S2,S3,U1,U3,W1,W2,W3,EHR[1:N],R,RF[1:5,1:N],
            ERR[1:3],A1,A2[1:N,1:N];
    integer array PS1,PS2[1:N];
    boolean START,TWO,HALV;

    procedure DFDY(T,Y,A); real T; array Y,A;
    begin integer I,J; real SL; array F1,F2[1:N];
        DERIV(T,Y,F1,N);
        for I:=1 step 1 until N do 
        begin 
            SL:="-6*Y[I]; if ABS(SL)<"-6 then SL:="-6;
            Y[I]:=Y[I]+SL; DERIV(T,Y,F2,N);
            for J:=1 step 1 until N do 
            A[J,I]:=(F2[J]-F1[J])/SL;
            Y[I]:=Y[I]-SL;
        end 
    end DFDY;

    procedure STARTV(Y,T); value T; real T; array Y;
    begin real A,B,C;
        A:=(T-T1)/(T1-T2); B:=(T-T2)/(T1-T3);
        C:=(T-T1)/(T2-T3)*B; B:=A*B;
        A:=1+A+B; B:=A+C-1;
        MULVEC(1,N,0,Y,S1,A); ELMVEC(1,N,0,Y,S2,-B);
        ELMVEC(1,N,0,Y,S3,C)
    end STARTV

    procedure ITERATE(Z,Y,A,H,T,WEIGHTS,FAIL,PS);
    array Z,Y,A,WEIGHTS; real H,T; label FAIL;
    integer array PS;
    begin integer IT,LIT; real MAX,MAX1,CONV; array DZ,F1[1:N];
        for I:=1 step 1 until N do Z[I]:=(Z[I]+Y[I])/2;
        IT:=LIT:=1; CONV:=1;
    ATER: DERIV(T,Z,F1,N);
        for I:=1 step 1 until N do 
        F1[I]:=DZ[I]:=Z[I]-H*F1[I]/2-Y[I];
        SOL(A,N,PS,DZ);
        ELMVEC(1,N,0,Z,DZ,-1);
        MAX:=0;
        for I:=1 step 1 until N do 
        MAX:=MAX+(WEIGHTS[I]*DZ[I])**2;
        MAX:=SQRT(MAX);
        if MAX*CONV<EPS/10 then goto OUT;
        IT:=IT+1; if IT=2 then goto ASS;
        CONV:=MAX/MAX1;
        if CONV>.2 then 
        begin if LIT=0 then goto FAIL;
            LIT:=0; CONV:=1; IT:=1;
            RECOMP(A,H,T,Z,FAIL,PS);
        end;
    ASS: MAX1:=MAX;
        goto ATER;
    OUT: for I:=1 step 1 until N do Z[I]:=2*Z[I]-Y[I];
    end ITERATE;

    procedure RECOMP(A,H,T,Y,FAIL,PS);
    real H,T; array A,Y; label FAIL; integer array PS;
    begin real SL; array AUX[1:3];
        SL:=H/2;
        if not AVAILABLE(T,Y,A,N) then DFDY(T,Y,A);
        for I:=1 step 1 until N do 
        begin MULROW(1,N,I,I,A,A,-SL); A[I,I]:=1+A[I,I]
        end;
        AUX[2]:="-14;
        DEC(A,N,AUX,PS);
        if AUX[3]<N then goto FAIL
    end RECOMP;

    procedure INITIALIZATION;
    begin H2:=HNEW; H:=H2/2;
        DUPVEC(1,N,0,S1,Y0); DUPVEC(1,N,0,S2,Y0); DUPVEC(1,N,0,S3,Y0);
        DUPVEC(1,N,0,W1,Y0); DUPROWVEC(1,N,1,R,Y0);
        INIVEC(1,N,U1,0); INIVEC(1,N,W2,0);
        INIMAT(2,5,1,N,R,0); INIMAT(1,5,1,N,RF,0);
        T:=T1:=T0; T2:=T0-2*H-"6; T3:=2*T2+1;
        RECOMP(A1,H,T,S1,MISS,PS1);RECOMP(A2,H2,T,W1,MISS,PS2);
    end 

    procedure ONE LARGE STEP;
    begin STARTV(Z,T+H);
        ITERATE(Z,S1,A1,H,T+H/2,WEIGHTS,MISS,PS1);
        DUPVEC(1,N,0,Y,Z);
        STARTV(Z,T+H2);
        ITERATE(Z,Y,A1,H,T+3*H/2,WEIGHTS,MISS,PS1);
        DUPVEC(1,N,0,U3,U1); DUPVEC(1,N,0,U1,Y);
        DUPVEC(1,N,0,S3,S2); DUPVEC(1,N,0,S2,S1);
        DUPVEC(1,N,0,S1,Z);
        ELMVEC(1,N,0,Z,W1,1); ELMVEC(1,N,0,Z,S2,-1);
        ITERATE(Z,W1,A2,H2,T+H,WEIGHTS,MISS,PS2);
        T3:=T2; T2:=T1; T1:=T+H2;
        DUPVEC(1,N,0,W3,W2); DUPVEC(1,N,0,W2,W1); DUPVEC(1,N,0,W1,Z);
    end;

    procedure CHANGE OF INFORMATION;
    begin real ALF1,C1,C2,C3; array KOF[2:4,2:4],E,D[1:4];
        integer I, K;
        C1:=HNEW/H2; C2:=C1*C1; C3:=C2*C1;
        KOF[2,2]:=C1; KOF[2,3]:=(C1-C2)/2; KOF[2,4]:=C3/6-C2/2+C1/3;
        KOF[3,3]:=C2; KOF[3,4]:=C2-C3; KOF[4,4]:=C3;
        for I:=1 step 1 until N do 
        U1[I]:=R[2,I]+R[3,I]/2+R[4,I]/3;
        ALF1:=MATVEC(1,N,1,RF,U1)/VECVEC(1,N,0,U1,U1);
        ALF:=(ALF+ALF1)*C1;
        for I:=1 step 1 until N do 
        begin 
            E[1]:=RF[1,I]-ALF1*U1[I];
            E[2]:=RF[2,I]-ALF1*2*R[3,I];
            E[3]:=RF[3,I]-ALF1*4*R[4,I];
            E[4]:=RF[4,I];
            D[1]:=R[1,I]; RF[1,I]:=E[1]:=E[1]*C2;
            for K:=2 step 1 until 4 do 
            begin R[K,I]:=D[K]:=MATMAT(K,4,K,I,KOF,R);
                RF[K,I]:=E[K]:=C2*MATVEC(K,4,K,KOF,E)
            end K;
            S1[I]:=D[1]+E[1];W1[I]:=D[1]+4*E[1];
            S2[I]:=S1[I]-(D[2]+E[2]/2);
            S3[I]:=S2[I]-(D[2]+E[2])+(D[3]+E[3]/2);
        end I;
        T3:=T-HNEW; T2:=T-HNEW/2; T1:=T;
        H2:=HNEW; H:=H2/2; ERR[1]:=0;
        if HALV then 
        begin for I:=1 step 1 until N do PS2[I]:= PS1[I];
            DUPMAT(1,N,1,N,A2,A1) end;
        if TWO then 
        begin for I:=1 step 1 until N do PS1[I]:= PS2[I];
            DUPMAT(1,N,1,N,A1,A2)
        end else RECOMP(A1,HNEW/2,T,S1,MISS,PS1);
        if ^HALV then RECOMP(A2,HNEW,T,W1,MISS,PS2);
    end HNEW^=H2
    procedure BACKWARD DIFFERENCES;
    for I:=1 step 1 until N do 
    begin real B0,B1,B2,B3;
        B1:=(U1[I]+2*S2[I]+U3[I])/4;
        B2:=(W1[I]+2*W2[I]+W3[I])/4;
        B3:=(S3[I]+2*U3[I]+S2[I])/4;
        B2:=(B2-B1)/3; B0:=B1-B2;
        B2:=B2-(S1[I]-2*S2[I]+S3[I])/16;
        B1:=2*B3-(B2+RF[1,I])-(B0+R[1,I])/2;
        B3:=0;
        for K:=1 step 1 until 4 do 
        begin B1:=B1-B3; B3:=R[K,I]; R[K,I]:=B0; B0:=B0-B1
        end; R[5,I]:=B0;
        for K:=1 step 1 until 4 do 
        begin B3:=RF[K,I]; RF[K,I]:=B2; B2:=B2-B3 end;
        RF[5,I]:=B2;
    end;

    procedure ERROR ESTIMATES;
    begin real C0,C1,C2,C3,B0,B1,B2,B3,W,SL1,SN,LR;
        C0:=C1:=C2:=C3:=0;
        for I:=1 step 1 until N do 
        begin W:=WEIGHTS[I]**2;
            B0:=RF[4,I]/36; C0:=C0+B0*B0*W; LR:=ABS(B0);
            B1:=RF[1,I]+ALF*R[2,I]; C1:=C1+B1*B1*W;
            B2:=RF[3,I]; C2:=C2+B2*B2*W;
            SL1:=ABS(RF[1,I]-RF[2,I]);
            SN:=if SL1<"-10 then 1else ABS(RF[1,I]-R[4,I]/6)/SL1;
            if SN>1 then SN:=1;
            if START then begin SN:=SN**4; LR:=LR*4 end;
            EHR[I]:=B3:=SN*EHR[I]+LR; C3:=C3+B3*B3*W;
        end I;
        B0:=ERR[1];
        ERR[1]:=B1:=SQRT(C0); ERR[2]:=SQRT(C1);
        ERR[3]:=SQRT(C3)+SQRT(C2)/2;
        LQ:=EPS/(if B0<B1 then B1"ELSE" B0);
        if B0<B1 and LQ>=80 then LQ:=10;
    end;

    procedure REJECT;
    if START then 
    begin HNEW:=LQ**(1/5)*H/2; goto INIT
    end else 
    begin for K:=1,2,3,4,1,2,3 do ELMROW(1,N,K,K+1,R,R,-1);
        for K:=1,2,3,4 do ELMROW(1,N,K,K+1,RF,RF,-1);
        T:=T-H2; HALV:=true; HNEW:=H; goto MSTP
    end 

    procedure STEPSIZE;
    if LQ<2 then 
    begin HALV:=true; HNEW:=H end else 
    begin if LQ>80 then 
        HNEW:=(if LQ>5120 then (LQ/5)**(1/5) else 2)*H2;
        if HNEW>HMAX then HNEW:=HMAX;
        if TEND>T and TEND-T<HNEW then HNEW:=TEND-T;
        TWO:=HNEW=2*H2;
    end;

    if PRESCH then H:=H0 else 
    begin if H0>HMAX then H:=HMAX else H:=H0;
        if H>(TEND-T0)/4 then H:=(TEND-T0)/4;
    end;
    HNEW:=H;
    ALF:=0; T:=TP:=T0;
    INIVEC(1,3,ERR,0); INIVEC(1,N,EHR,0);
    DUPROWVEC(1,N,1,R,Y0); INIMAT(2, 5, 1, N, R, 0.0);
    CONTROL(TP,T,H,HNEW,R,ERR,N);
INIT: INITIALIZATION; START:=true;
    for ECI:=0,1,2,3 do 
    begin ONE LARGE STEP; T:=T+H2;
        if ECI>0 then 
        begin BACKWARD DIFFERENCES; UPDATE(WEIGHTS,S2,N) end 
    end;
    ECI:=4;
MSTP: if HNEW^=H2 then 
    begin ECI:=1; CHANGE OF INFORMATION;
        ONE LARGE STEP; T:=T+H2; ECI:=2;
    end;
    ONE LARGE STEP;
    BACKWARD DIFFERENCES;
    UPDATE(WEIGHTS,S2,N);
    ERROR ESTIMATES;
    if ECI<4 and LQ>80 then LQ:=20;
    HALV:=TWO:=false;
    if PRESCH then goto TRYCK;
    if LQ<1 then REJECT else STEPSIZE;
TRYCK: if TP<=T then CONTROL(TP,T,H,HNEW,R,ERR,N);
    if START then START:=false;
    if HNEW=H2 then T:=T+H2; ECI:=ECI+1;
    if T<TEND+H2 then goto MSTP else goto END;
MISS: FAIL:=PRESCH;
    if ^ FAIL then 
    begin if ECI>1 then T:=T-H2;
        HALV:=TWO:=false; HNEW:=H2/2;
        if START then goto INIT else goto TRYCK
    end;
END:
end IMPEX;
        eop