code 33017 ;
procedure RK4NA(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, IV, IV0;
boolean FIR, FIRST, REJ;
real H, COND0, COND1, FHM, ABSH, TOL, FH, MAX, X0,
X1, S, HMIN, HL, MU, MU1;
array XL, DISCR, Y[0:N], K[0:5,0:N], E1[1:2];
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:= 1 step 1 until N do Y[J]:= FXJ;
P:= H / Y[IV];
for I:= 0 step 1 until N do
begin if I ^= IV then K[T,I]:= Y[I] * P end
end F;
if D = 2 then goto INTEGRATE; if D = 3 then
begin for I:= 0 step 1 until N do X[I]:= XL[I];
F(0)
end
else if D = 1 then
begin real P;
P:= H / Y[IV];
for I:= 0 step 1 until N do
begin if I ^= IV then K[0,I]:= P * Y[I] end
end
else
for I:= 0 step 1 until N do
begin if I ^= IV then K[0,I]:= K[0,I] * MU end;
for I:= 0 step 1 until N do X[I]:= XL[I] + (if I
= IV then H else K[0,I]) / 4.5; F(1);
for I:= 0 step 1 until N do X[I]:= XL[I] + (if I
= IV then H * 4 else (K[0,I] + K[1,I] * 3)) / 12;
F(2);
for I:= 0 step 1 until N do X[I]:= XL[I] + (if I
= IV then H * .5 else (K[0,I] + K[2,I] * 3) / 8);
for I:= 0 step 1 until N do X[I]:= XL[I] + (if I
= IV then H * .8 else (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] +
(if I = IV then H else (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
begin if I ^= IV then DISCR[I]:= ABS(K[0,I] * 21
- K[2,I] * 162 + K[3,I] * 224 - K[4,I] *
125 + K[5,I] * 42) / 14
end;
goto END
end;
INTEGRATE: for I:= 0 step 1 until N do X[I]:= XL[I]
+ (if I = IV then H else ( - 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
begin if I ^= IV then X[I]:= XL[I] + (K[0,I] * 35
+ K[2,I] * 162 + K[4,I] * 125 + K[5,I] * 14) / 336
end ;
END ..
end RKSTEP ;
real procedure FZERO;
begin if S = X0 then FZERO:= COND0 else if S = X1
then FZERO:= COND1 else
begin RKSTEP(S - XL[IV], 3); FZERO:= B end
end FZERO;
if FI then
begin for I:= 0 step 1 until N do D[I + 3]:= XA[I];
D[0]:= D[2]:= 0
end;
D[1]:= 0;
for I:= 0 step 1 until N do X[I]:= XL[I]:= D[I + 3];
IV:= D[0]; H:= D[2]; FIRST:= FIR:= true; Y[0]:= 1;
goto CHANGE;
AGAIN: ABSH:= ABS(H); if ABSH < HMIN then
begin H:= if H > 0 then HMIN else - HMIN;
ABSH:= ABS(H)
end;
RKSTEP(H, I); REJ:= false; FHM:= 0;
for I:= 0 step 1 until N do
begin if I ^= IV then
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 if ABSH <= HMIN then
begin for I:= 0 step 1 until N do
begin if I ^= IV then X[I]:= XL[I] + K[0,I]
else X[I]:= XL[I] + H
end;
D[1]:= D[1] + 1; FIRST:= true; goto NEXT
end;
H:= H * MU; I:= 0; goto AGAIN
end;
if FIRST then
begin FIRST:= FIR; HL:= H; H:= MU * H; goto ACCEPT
end;
FH:= MU * H / HL + MU - MU1; HL:= H; H:= FH * H;
ACCEPT: RKSTEP(HL, 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;
for I:= 0 step 1 until N do D[I + 3]:= XL[I]:= X[I];
CHANGE: IV0:= IV;
for J:= 1 step 1 until N do Y[J]:= FXJ;
MAX:= ABS(Y[IV]);
for I:= 0 step 1 until N do
begin if ABS(Y[I]) > MAX then
begin MAX:= ABS(Y[I]); IV:= I end
end;
if IV0 ^= IV then
begin FIRST:= true; D[0]:= IV;
D[2]:= H:= Y[IV] / Y[IV0] * H
end;
X0:= XL[IV]; if FIR then
begin HMIN:= E[0] + E[1];
for I:= 1 step 1 until N do
begin H:= E[2 * I] + E[2 * I + 1];
if H < HMIN then HMIN:= H
end;
H:= E[2 * IV] + E[2 * IV + 1];
if (FI and (Y[L]/Y[IV]*H<0 eqv POS)) or
(not FI and D[2]*H<0) then H:= -H
end;
I:= 1; goto AGAIN;
ZERO: E1[1]:= E[2 * N + 2]; E1[2]:= E[2 * N + 3];
X1:=X[IV] ; S:=X0 ;
ZEROIN(S,X1,FZERO,ABS(E1[1]*S) + ABS(E1[2])) ; X0:=S ; X1:=X[IV];
RKSTEP(X0 - XL[IV], 3);
for I:= 0 step 1 until N do D[I + 3]:= X[I]
end RK4NA;
eop