code 33018 ;
   procedure RK5NA(X, XA, B, FXJ, J, E, D, FI, N, L, POS);
   value FI, N, L, POS; integer J, N, L; boolean FI, POS;
   real B, FXJ; array X, XA, E, D;
   begin integer I;
      boolean FIRST, FIR, REJ;
      real FHM, S, S0, COND0, S1, COND1, H, ABSH, TOL, FH,
      HL, MU, MU1;
      array Y, XL, DISCR[0:N], K[0:5,0:N], E1[1:2];
      real procedure SUM(J,A,B,XJ) ; integer J,A,B ; real XJ ;
      begin real S ; S:= 0 ;
         for J:=A step 1 until B do S:=S+XJ ; SUM:= S
      end SUM ;
      procedure RKSTEP(H, D); value H, D; integer D; real H;
      begin integer I;
         procedure F(T); value T; integer T;
         begin integer I;
            real P;
            for J:= 0 step 1 until N do Y[J]:= FXJ;
            P:= H / SQRT(SUM(I, 0, N, Y[I] ** 2));
            for I:= 0 step 1 until N do K[T,I]:= Y[I] * P
         end F;
         if D = 2 then goto INTEGRATE; if D = 1 then 
         begin for I:= 0 step 1 until N do K[0,I]:= K[0,I]
            * MU; goto A
         end;
         for I:= 0 step 1 until N do X[I]:= XL[I]; F(0);
      A: for I:= 0 step 1 until N do X[I]:= XL[I] +
         K[0,I] / 4.5; F(1);
         for I:= 0 step 1 until N do X[I]:= XL[I] + (K[0,I]
         + K[1,I] * 3) / 12; F(2);
         for I:= 0 step 1 until N do X[I]:= XL[I] + (K[0,I]
         + K[2,I] * 3) / 8; F(3);
         for I:= 0 step 1 until N do X[I]:= XL[I] + (K[0,I]
         * 53 - K[1,I] * 135 + K[2,I] * 126 + K[3,I] * 56)
         / 125; F(4); if D <= 1 then 
         begin for I:= 0 step 1 until N do X[I]:= XL[I] +
            (K[0,I] * 133 - K[1,I] * 378 + K[2,I] * 276 +
            K[3,I] * 112 + K[4,I] * 25) / 168; F(5);
            for I:= 0 step 1 until N do DISCR[I]:=
            ABS(K[0,I] * 21 - K[2,I] * 162 + K[3,I] * 224
            - K[4,I] * 125 + K[5,I] * 42) / 14; goto END
         end;
      INTEGRATE: for I:= 0 step 1 until N do X[I]:= XL[I]
         + ( - K[0,I] * 63 + K[1,I] * 189 - K[2,I] * 36 -
         K[3,I] * 112 + K[4,I] * 50) / 28; F(5);
         for I:= 0 step 1 until N do X[I]:= XL[I] + (K[0,I]
         * 35 + K[2,I] * 162 + K[4,I] * 125 + K[5,I] * 14)
         / 336;
      END:
      real procedure FZERO;
      begin if S = S0 then FZERO:= COND0 else if S = S1
         then FZERO:= COND1 else 
         begin RKSTEP(S - S0, 3); FZERO:= B end 
      end FZERO;

      if FI then 
      begin for I:= 0 step 1 until N do D[I + 3]:= XA[I];
         D[1]:= D[2]:= 0
      end;
      for I:= 0 step 1 until N do X[I]:= XL[I]:= D[I + 3];
      S:= D[1]; FIRST:= FIR:= true; H:= E[0] + E[1];
      for I:= 1 step 1 until N do 
      begin ABSH:= E[2 * I] + E[2 * I + 1];
         if H > ABSH then H:= ABSH
      end;
      if FI then 
      begin J:= L; if FXJ * H < 0 eqv POS then H:= - H end 
      else if D[2] * H < 0 then H:= - H; I:= 0;
   AGAIN: RKSTEP(H, I); REJ:= false; FHM:= 0;
      ABSH:= ABS(H);
      for I:= 0 step 1 until N do 
      begin TOL:= E[2 * I] * ABS(K[0,I]) + E[2 * I + 1] *
         ABSH; REJ:= TOL < DISCR[I] or REJ;
         FH:= DISCR[I] / TOL; if FH > FHM then FHM:= FH
      end;
      MU:= 1 / (1 + FHM) + .45; if REJ then 
      begin H:= H * MU; I:= 1; goto AGAIN end;
      if FIRST then 
      begin FIRST:= FIR; HL:= H; H:= MU * H end 
      else 
      begin FH:= MU * H / HL + MU - MU1; HL:= H; H:= FH * H
      end;
   ACCEPT: RKSTEP(HL, 2); MU1:= MU; S:= S + HL;
      if FIR then 
      begin COND0:= B; FIR:= false; if not FI then H:= D[2]
      end 
      else 
      begin D[2]:= H; COND1:= B;
         if COND0 * COND1 <= 0 then goto ZERO;
         COND0:= COND1
      end;
      for I:= 0 step 1 until N do D[I + 3]:= XL[I]:= X[I];
      D[1]:= S0:= S; I:= 0; goto AGAIN;
   ZERO: E1[1]:= E[2 * N + 2]; E1[2]:= E[2 * N + 3];
      S1:=S ; S:=S0 ;
      ZEROIN(S,S1,FZERO,ABS(E1[1]*S)+ABS(E1[2])) ;
      RKSTEP(S - S0, 3);
      for I:= 0 step 1 until N do D[I + 3]:= X[I]; D[1]:= S
   end RK5NA;
        eop