code 33016 ;
   procedure RK4A(X, XA, B, Y, YA, FXY, E, D, FI, XDIR,
   POS); value FI, XDIR, POS; boolean FI, XDIR, POS;
   real X, XA, B, Y, YA, FXY; array E, D;
   begin integer I;
      boolean IV, FIRST, FIR, REJ;
      real K0, K1, K2, K3, K4, K5, FHM, ABSH, DISCR, S, XL,
      COND0, S1, COND1, YL, HMIN, H, ZL, TOL, HL, MU, MU1;
      array E1[1:2];

      procedure RKSTEP(X, XL, H, Y, YL, ZL, FXY, D);
      value XL, YL, ZL, H; real X, XL, H, Y, YL, ZL, FXY;
      integer D;
      begin if D = 2 then goto INTEGRATE; if D = 3 then 
         begin X:= XL; Y:= YL; K0:= FXY * H end 
         else if D = 1 then K0:= ZL * H else K0:= K0 * MU;
         X:= XL + H / 4.5; Y:= YL + K0 / 4.5; K1:= FXY * H;
         X:= XL + H / 3; Y:= YL + (K0 + K1 * 3) / 12;
         K2:= FXY * H; X:= XL + H * .5;
         Y:= YL + (K0 + K2 * 3) / 8; K3:= H * FXY;
         X:= XL + H * .8;
         Y:= YL + (K0 * 53 - K1 * 135 + K2 * 126 + K3 *
         56) / 125; K4:= FXY * H; if D <= 1 then 
         begin X:= XL + H;
            Y:= YL + (K0 * 133 - K1 * 378 + K2 * 276 + K3
            * 112 + K4 * 25) / 168; K5:= FXY * H;
            DISCR:= ABS(K0 * 21 - K2 * 162 + K3 * 224 - K4
            * 125 + K5 * 42) / 14; goto END
         end;
      INTEGRATE: X:= XL + H;
         Y:= YL + ( - K0 * 63 + K1 * 189 - K2 * 36 - K3 *
         112 + K4 * 50) / 28; K5:= FXY * H;
         Y:= YL + (K0 * 35 + K2 * 162 + K4 * 125 + K5 *
         14) / 336;
      END:
      real procedure FZERO;
      begin if IV then 
         begin if S = XL then FZERO:= COND0 else if S = S1
            then FZERO:= COND1 else 
            begin RKSTEP(X, XL, S - XL, Y, YL, ZL, FXY, 3);
               FZERO:= B
            end 
         end 
         else 
         begin if S = YL then FZERO:= COND0 else if S = S1
            then FZERO:= COND1 else 
            begin RKSTEP(Y, YL, S - YL, X, XL, ZL, 1 /
               FXY, 3); FZERO:= B
            end 
         end 
      end FZERO;

      if FI then 
      begin D[3]:= XA; D[4]:= YA; D[0]:= 1 end;
      D[1]:= 0; X:= XL:= D[3]; Y:= YL:= D[4]; IV:= D[0] > 0;
      FIRST:= FIR:= true; HMIN:= E[0] + E[1];
      H:= E[2] + E[3]; if H < HMIN then HMIN:= H;
   CHANGE: ZL:= FXY; if ABS(ZL) <= 1 then 
      begin if not IV then 
         begin D[2]:= H:= H / ZL; D[0]:= 1;
            IV:= FIRST:= true 
         end;
         if FIR then goto A; I:= 1; goto AGAIN
      end 
      else 
      begin if IV then 
         begin if not FIR then D[2]:= H:= H * ZL; D[0]:= - 1;
            IV:= false; FIRST:= true 
         end;
         if FIR then 
         begin H:= E[0] + E[1];
         A: if (if FI then (if IV eqv XDIR then H else 
            H * ZL) < 0 eqv POS else H * D[2] < 0) then H:= - H
         end;
         I:= 1
   AGAIN: ABSH:= ABS(H); if ABSH < HMIN then 
      begin H:= SIGN(H) * HMIN; ABSH:= HMIN end;
      if IV then 
      begin RKSTEP(X, XL, H, Y, YL, ZL, FXY, I);
         TOL:= E[2] * ABS(K0) + E[3] * ABSH
      end 
      else 
      begin RKSTEP(Y, YL, H, X, XL, 1 / ZL, 1 / FXY, I);
         TOL:= E[0] * ABS(K0) + E[1] * ABSH
      end;
      REJ:= DISCR > TOL; MU:= TOL / (TOL + DISCR) + .45;
      if REJ then 
      begin if ABSH <= HMIN then 
         begin if IV then 
            begin X:= XL + H; Y:= YL + K0 end 
            else 
            begin X:= XL + K0; Y:= YL + H end;
            D[1]:= D[1] + 1; FIRST:= true; goto NEXT
         end;
         H:= H * MU; I:= 0; goto AGAIN
      end REJ;
      if FIRST then 
      begin FIRST:= FIR; HL:= H; H:= MU * H; goto ACCEPT
      end;
      FHM:= MU * H / HL + MU - MU1; HL:= H; H:= FHM * H;
   ACCEPT: if IV then RKSTEP(X, XL, HL, Y, YL, ZL, FXY,
      2) else RKSTEP(Y, YL, HL, X, XL, ZL, 1 / FXY, 2);
      MU1:= MU;
   NEXT: if FIR then 
      begin FIR:= false; COND0:= B;
         if not (FI or REJ) then H:= D[2]
      end 
      else 
      begin D[2]:= H; COND1:= B;
         if COND0 * COND1 <= 0 then goto ZERO;
         COND0:= COND1
      end;
      D[3]:= XL:= X; D[4]:= YL:= Y; goto CHANGE;
   ZERO: E1[1]:= E[4]; E1[2]:= E[5];
      S1:= if IV then X else Y;
      S:= if IV then XL else YL ;
      ZEROIN(S,S1,FZERO,ABS(E1[1]*S)+ABS(E1[2])) ;
      S1:= if IV then X else Y ;
      if IV then RKSTEP(X, XL, S - XL, Y, YL, ZL, FXY, 3)
      else RKSTEP(Y, YL, S - YL, X, XL, ZL, 1 / FXY,
      3); D[3]:= X; D[4]:= Y
   end RK4A;
        eop