{   26/01/84 011141605  HLL     RKF45    PASCAL             PASCAL   }

const
   NEQN = 20;

type
   workarray = array[1..NEQN] of real;  

procedure rkf45

   (   n       : integer;   { number of equations }
    function
       f                    { compute y_prime = dy / dx;
                              return true if ok; else false }
        (   n  : integer;   { number of equations }
            x  : real;      { independent variable }
            y  : workarray; { dependent variable }
         var
            y_prime : workarray  { dy / dx }
        )
               : boolean;
    var
       y_prime : workarray; { derivatives of y computed by f }
       x0      : real;      { initial value of independent variable }
    var
       y0      : workarray; { initial value of dependent variable }
       x1      : real;      { final value of independent variable }
       epsilon : workarray; { integration tolerance, expressed as an
                              absolute value, not a fraction }
       firsth,              { initial step size }
       minh                 { smallest step size tolerated }
               : real;  
    var
       y1      : workarray; { final value of dependent variable; only
                              returned if success is true }
    var
       lasth   : real;      { last value of step size; recommended for
                              restart if desired }
    var
       failure : integer    { was integration successful? }
   ); 

    {-----------------------------------------------------------------
    |  Runge-Kutta-Fehlberg integration of the system of n ordinary   |
    |  differential equations                                         |
    |                                                                 |
    |    dy(x) [i]                                                    |
    |    -----      = y_prime(x) [i];                                 |
    |    dx                                                           |
    |                                                                 |
    |    y(x0) [i]  = y0 [i]; --> y(x1) [i] = y1 [i];                 |
    |                                                                 |
    |            (i = 1...n)                                          |
    |                                                                 |
    |  using the 5-stage fourth order Runge-Kutta method of Fehlberg. |
    |                                                                 |
    |  The function f(n, x, y, y_prime) computes y_prime from the     |
    |  values of x and y(x) and returns true if the calculation is    |
    |  successful, else false.                                        |
    |                                                                 |
    |  In application rkf45 should be used in a series of steps       |
    |  to cover the desired total interval; the method will take      |
    |  finer steps within this interval as required.  Internally      |
    |  the step size is halved if the error estimate per unit step    |
    |  exceeds epsilon.   The step size is doubled if the error       |
    |  estimate is less than epsilon / 32, unless this exceeds the    |
    |  range available.  The process terminates if the step size      |
    |  becomes less than minh.  The variable "failure" is 0 if the    |
    |  integration is successful, 1 if the number of equations        |
    |  exceeds NEQN, 2 if the step size became less than minh without |
    |  reaching the requested error tolerance, 3 if the derivative    |
    |  calculation failed.                                            |
    |                                                                 |
    |  Taken from An Introduction to Numerical Methods with Pascal    |
    |  p. 283 by L.V. Atkinson and P.J. Harley, (Addison Wesley).     |
    |                                                                 |
    |  Modified                                                       |
    |  (1) to protect the routine from bad step size and collapsed    |
    |      x interval,                                                |
    |  (2) to allow stepping in the  negative direction,              |
    |  (3) to return several error status indications,                |
    |  (4) to work in n-dimensions.                                   |
    |  (5) better value of nearzero (roundoff handling)               |
    |  (6) pass y0 and y_prime as var parameters.                     |
    |                                                                 |
    |                     26 January 1984                             |
    |                        H. L. Lynch                              |
    |                                                                 |
    ------------------------------------------------------------------}

   const
      mach_acc = 1.0e-15;       { relative machine precision for real }
   var
      yt,                       { temporary storage for calc. kn }
      k1,       k2,             { interpolated devivative values }
      k3,       k4,             { interpolated derivative values }
      k5,       k6              { interpolated derivative values }
         : workarray;
      nearzero,                 { within machine accuracy of end point }
      error,                    { error estimate }
      x,                        { curr. value of independent variable}
      h,                        { current step size }
      ah,                       { absolute value of current step size}
      tolerance                 { temp. var. in error tolerance calc.}
         : real;  
      state                     { current state of program }
         : (stepping, x1reached, htoosmall, bad_derivative);
      good_deriv,               { good derivative returned from f() }
      err_too_large,            { need to reduce   step size }
      err_too_small,            { need to increase step size }
      shortrange                { flag to prevent step beyond x1 and
                                  return good estimated h }
         : boolean;
      i                         { dummy index for arrays }
         : integer;
begin      { Runge-Kutta-Fehlberg rkf45 }

   if (n > NEQN) then    { Attempt to over run internal arrays }
   begin
      failure := 1;
      return;
   end;

   x := x0;
   for i := 1 to n do
      y1[i] := y0[i];
   shortrange := false;
   state := stepping;
   h := firsth;

   { check for bad input }

   if ((h * (x1 - x0)) < 0.0) then  { step is in wrong direction }
      h := -h;
   if ((((x1 - x0) / h) < 1) or (h = 0)) then   { bad step size }
   begin
      h := x1 - x0;
      if (abs(h) < minh) then     { interval has collapsed }
         return;
   end;

   ah := abs(h);
   nearzero := mach_acc * abs(x1 - x0);

   repeat { while state = stepping, checked at the end }
      good_deriv := f(n, x, y1, y_prime);
      for i := 1 to n do
         k1[i] := y_prime[i];

      for i := 1 to n do
         yt[i] := y1[i] + h * k1[i] / 4.0;
      good_deriv := f(n, x + h / 4.0, yt, y_prime) and good_deriv;
      for i := 1 to n do
         k2[i] := y_prime[i];

      for i := 1 to n do
         yt[i] := y1[i] + h * (3.0 * k1[i] + 9.0 * k2[i]) / 32.0;
      good_deriv := f(n, x + 3.0 * h / 8.0, yt, y_prime)
                     and good_deriv;
      for i := 1 to n do
         k3[i] := y_prime[i];

      for i := 1 to n do
         yt[i] := y1[i] + h * (1932.0 * k1[i]
                  - 7200.0 * k2[i] + 7296.0 * k3[i]) / 2197.0;
      good_deriv := f(n, x + 12.0 * h / 13.0, yt, y_prime)
                     and good_deriv;
      for i := 1 to n do
         k4[i] := y_prime[i];

      for i := 1 to n do
         yt[i] := y1[i] + h * (439.0 * k1[i] / 216.0 - 8.0 * k2[i]
                  + 3680.0 * k3[i] / 513.0 - 845.0 * k4[i] / 4104.0);
      good_deriv := f(n, x + h, yt, y_prime) and good_deriv;
      for i := 1 to n do
         k5[i] := y_prime[i];

      for i := 1 to n do
         yt[i] := y1[i] + h * (-8.0 * k1[i] / 27.0 + 2.0 * k2[i]
                  - 3544.0 * k3[i] / 2565.0 + 1859.0 * k4[i] / 4104.0
                  - 11.0 * k5[i] / 40.0);
      good_deriv := f(n, x + h / 2.0, yt, y_prime) and good_deriv;
      for i := 1 to n do
         k6[i] := y_prime[i];

      if (not good_deriv) then
         state := bad_derivative;

      { Set up err_too_large := at least one error is too large;
               err_too_small := all errors are too small.
      }

      err_too_large := false;
      err_too_small := true;
      for i := 1 to n do
      begin
         error := ah * abs(k1[i] / 360.0 - 128.0 * k3[i] / 4275.0
                           - 2197.0 * k4[i] / 75240.0
                           + k5[i] / 50.0 + 2 * k6[i] / 55.0);
         tolerance := ah * epsilon[i];
         err_too_large := (error > tolerance)          or err_too_large;
         err_too_small := (error <= tolerance / 32.0) and err_too_small;
      end; { end i loop }

      if (err_too_large) then
      begin
         h := h / 2.0;
         ah := abs(h);
         if (ah < minh) then
            state := htoosmall;
      end
      else
      begin  { Good result, take the step. }
         x := x + h;
         for i := 1 to n do
         begin
            y1[i] := y1[i] + h * (25.0 * k1[i] / 216.0
                     + 1408.0 * k3[i] / 2565.0
                     + 2197.0 * k4[i] / 4104.0 - k5[i] / 5.0);
         end;
         if (err_too_small) then
         begin
            h := h * 2.0;
            ah := abs(h);
         end;
         if (abs(x1 - x) <= nearzero) then
            state := x1reached
         else
         begin
            shortrange := (ah > abs(x1 - x));
            if (shortrange) then
            begin
               lasth := h;
               h := x1 - x;
               ah := abs(h);
            end;
         end;  { end x1 not yet reached }
      end;  { end err not too large }
   until (state <> stepping);

   if (not shortrange) then
      lasth := h;
   if      (state = x1reached)       then failure := 0
   else if (state = htoosmall)       then failure := 2
   else if (state = bad_derivative)  then failure := 3;
end { end Runge-Kutta-Fehlberg rkf45 } ;
