code 33180;
procedure DIFFSYS(X,XE,N,Y,DERIVATIVE,AETA,RETA,S,H0,OUTPUT);
value N;
integer N;
real X,XE,AETA,RETA,H0;
array Y,S;
procedure DERIVATIVE,OUTPUT;
begin real A,B,B1,C,G,H,U,V,TA,FC; integer I,J,K,KK,JJ,L,M,R,SR;
    array YA,YL,YM,DY,DZ[1:N],DT[1:N,0:6],D[0:6],YG,YH[0:7,1:N];
    boolean KONV,B0,BH,LAST;
    LAST:=false; H:=H0;
NEXT: if H*1.1>=XE-X then 
    begin LAST:=true; H0:=H; H:=XE-X+"-13 end;
    DERIVATIVE(X,Y,DZ); BH:=false;
    for I:=1 step 1 until N do YA[I]:=Y[I];
ANF: A:=H+X; FC:=1.5; B0:=false; M:=1; R:=2; SR:=3; JJ:=-1;
    for J:=0 step 1 until 9 do 
    begin if B0 then 
        begin D[1]:=16/9; D[3]:=64/9; D[5]:=256/9 end 
        else begin D[1]:=9/4; D[3]:=9; D[5]:=36 end;
        KONV:=true;
        if J>6 then begin L:=6; D[6]:=64; FC:=.6*FC end 
        else begin L:=J; D[L]:=M*M end;
        M:=M*2; G:=H/M; B:=G*2;
        if BH and J<8 then 
        begin for I:=1 step 1 until N do 
            begin YM[I]:=YH[J,I]; YL[I]:=YG[J,I] end 
        end 
        else 
                KK:=(M-2)/2; M:=M-1;
            for I:=1 step 1 until N do 
            begin YL[I]:=YA[I]; YM[I]:=YA[I]+G*DZ[I] end;
            for K:=1 step 1 until M do 
            begin DERIVATIVE(X+K*G,YM,DY);
                for I:=1 step 1 until N do 
                begin U:=YL[I]+B*DY[I]; YL[I]:=YM[I]; YM[I]:=U;
                    U:=ABS(U); if U>S[I] then S[I]:=U
                end;
                if K=KK and K^=2 then 
                begin JJ:=JJ+1; for I:=1 step 1 until N do 
                    begin YH[JJ,I]:=YM[I]; YG[JJ,I]:=YL[I] end 
                end 
            end 
        end;
        DERIVATIVE(A,YM,DY);
        for I:=1 step 1 until N do 
        begin V:=DT[I,0]; TA:=C:=DT[I,0]:=(YM[I]+YL[I]+G*DY[I])/2;
            for K:=1 step 1 until L do 
            begin B1:=D[K]*V; B:=B1-C; U:=V;
                if B^=0 then 
                begin B:=(C-V)/B; U:=C*B; C:=B1*B end;
                V:=DT[I,K]; DT[I,K]:=U; TA:=U+TA
            end;
            if ABS(Y[I]-TA)>RETA*S[I]+AETA then KONV:=false;
            Y[I]:=TA
        end;
        if KONV then goto END;
        D[2]:=4; D[4]:=16; B0:=^B0; M:=R; R:=SR; SR:=M*2
    end;
    BH:=^BH; LAST:=false; H:=H/2; goto ANF;
END: H:=FC*H; X:=A; OUTPUT; if not LAST then goto NEXT;
end DIFFSYS;
        eop