code 33033; procedure RKE (X, XE, N, Y, DER, DATA, FI, OUT); value FI, N; integer N; real X, XE; boolean FI; array Y, DATA; procedure DER, OUT; begin integer J; real XT, H, HMIN, INT, HL, HT, ABSH, FHM, DISCR, TOL, MU, MU1, FH, E1, E2; boolean LAST, FIRST, REJECT; array K0, K1, K2, K3, K4[1:N]; if FI then begin DATA[3]:= XE - X; DATA[4]:= DATA[5]:= DATA[6]:= 0 end; ABSH:= H:= ABS(DATA[3]); if XE < X then H:= - H; INT:= ABS(XE - X); HMIN:= INT * DATA[1] + DATA[2]; E1:= 12 * DATA[1] / INT; E2:= 12 * DATA[2] / INT; FIRST:= true; REJECT:= false; if FI then begin LAST:= true; goto STEP end; TEST: ABSH:= ABS(H); if ABSH < HMIN then begin H:= SIGN (XE - X) * HMIN; ABSH:= HMIN end; if H >= XE - X eqv H >= 0 then begin LAST:= true; H:= XE - X; ABSH:= ABS(H) end else LAST:= false; STEP: if ^ REJECT then begin for J:= 1 step 1 until N do K0[J]:= Y[J]; DER(X, K0) end; HT:= .184262134833347 * H; XT:= X + HT; for J:= 1 step 1 until N do K1[J]:= K0[J] * HT + Y[J]; DER(XT, K1); HT:= .690983005625053"-1 * H; XT:= 4 * HT + X; for J:= 1 step 1 until N do K2[J]:= (3 * K1[J] + K0[J]) * HT + Y[J]; DER(XT, K2); XT:= .5 * H + X; HT:= .1875 * H; for J:= 1 step 1 until N do K3[J]:=((1.74535599249993 * K2[J] - K1[J]) * 2.23606797749979 + K0[J]) * HT + Y[J]; DER(XT, K3); XT:= .723606797749979 * H + X; HT:= .4 * H; for J:= 1 step 1 until N do K4[J]:= (((.517595468166681 * K0[J] - K1[J]) * .927050983124840 + K2[J]) * 1.46352549156242 + K3[J]) * HT + Y[J]; DER(XT, K4); XT:= if LAST then XE else X + H; HT:= 2 * H; for J:= 1 step 1 until N do K1[J]:= ((((2 * K4[J] + K2[J]) * .412022659166595 + K1[J]) * 2.23606797749979 - K0[J]) * .375 - K3[J]) * HT + Y[J]; DER(XT, K1); REJECT:= false; FHM:= 0; for J:= 1 step 1 until N do begin DISCR:= ABS((1.6 * K3[J] - K2[J] - K4[J]) * 5 + K0[J] + K1[J]); TOL:= ABS(K0[J]) * E1 + E2; REJECT:= DISCR > TOL or REJECT; FH:= DISCR / TOL; if FH > FHM then FHM:= FH end; MU:= 1 / (1 + FHM) + .45; if REJECT then begin DATA[5]:= DATA[5] + 1; if ABSH <= HMIN then begin DATA[6]:= DATA[6] + 1; HL:= H; REJECT:= false; FIRST:= true; goto NEXT end; H:= MU * H; goto TEST end; if FIRST then begin FIRST:= false; HL:= H; H:= MU * H; goto ACC end; FH:= MU * H / HL + MU - MU1; HL:= H; H:= FH * H; ACC: MU1:= MU; HT:= HL / 12; for J:= 1 step 1 until N do Y[J]:= ((K2[J] + K4[J]) * 5 + K0[J] + K1[J]) * HT + Y[J]; NEXT: DATA[3]:= HL; DATA[4]:= DATA[4] + 1; X:= XT; OUT; if X ^= XE then goto TEST end RKE; eop