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;