procedure RK (x,y,n,FKT,eps,eta,xE,yE,fi); value x,y; integer n; Boolean fi; real x,eps,eta,xE; array y,yE; procedure FKT; comment RK integrates the system y’k=fk(x,y1,y2,...,yn)(k=1,2,...n) of differential equations with the method of Runge-Kutta with automatic search for appropriate length of integration step. Parameters are: The initial values x and y[k] for x and the unknown functions yk(x). The order n of the system. The procedure FKT(x,y,n,z) which represents the system to be integrated, i.e. the set of functions fk. The tolerance values eps and eta which govern the accuracy of the numerical integration. The end of the integration interval xE. The output parameter yE which represents the solution x=xE. The Boolean variable fi, which must always be given the value true for an isolated or first entry into RK. If however the functions y must be available at several meshpoints x0,x1,...,xn, then the procedure must be called repeatedly (with x=xk, xE=x(k+1), for k=0,1,...,n-1) and then the later calls may occur with fi=false which saves computing time. The input parameters of FKT must be x,y,z,n, the output parameter z represents the set of derivatives z[k]=fk(x,y[1],y[2],...,y[n]) for x and the actual y’s. A procedure comp enters as a non-local identifier; begin array z,y1,y2,y3[1:n]; real x1,x2,x3,H; Boolean out; integer k,j; own real s,Hs; procedure RK1ST (x,y,h,xe,ye); real x,h,xe; array y,ye; comment RK1ST integrates one single Runge-Kutta step with initial values x, y[k] which yields the output parameters xe=x+h and ye[k], the latter being the solution at xe. Important: the parameters n, FKT, z enter RK1ST as nonlocal entities; begin array w[1:n], a[1:5]; integer k,j; a[1] ≔ [2] ≔ [5] ≔ /2; a[3] ≔ [4] ≔ ; xe ≔ ; for k ≔ step 1 until n do ye[k] ≔ [k] ≔ [k]; for j ≔ step 1 until 4 do begin FKT(xe,w,n,z); xe ≔ +a[j]; for k ≔ step 1 until n do begin w[k] ≔ [k]+a[j]×z[k]; ye[k] ≔ e[k]+a[j+1]×z[k]/3 end k end j end RK1ST; Begin of program: if fi then begin H ≔ E-x; s ≔ end else H ≔ s; out ≔ alse; AA: if ((x+2·01×H-xE)>0) ≡ (H>0) then begin Hs ≔ ; out ≔ rue; H ≔ xE-x)/2 end if; RK1ST (x,y,2×H,x1,y1); BB: RK1ST (x,y,H,x2,y2); RK1ST (x2,y2,H,x3,y3); for k ≔ step 1 until n do if comp (y1[k],y3[k],eta) > eps then goto CC; comment comp(a,b,c) is a function designator, the value of which is the absolute value of the difference of the mantissae of a and b, after the exponents of these quantities have been made equal to the largest of the exponents of the originally given parameters a, b, c; x ≔ 3; if out then goto DD; for k ≔ step 1 until n do y[k] ≔ 3[k]; if s=5 then begin s ≔ ; H ≔ ×H end if; s ≔ +1; goto AA; CC: H ≔ ·5×X; out ≔ alse; x1 ≔ 2; for k ≔ step 1 until n do y1[k] ≔ 2[k]; goto BB; DD: for k ≔ step 1 until n do yE[k] ≔ 3[k] end RK;