{   11/01/84 011141605  HLL     TSTRKF45 PASCAL             PASCAL   }

program tstrkf45(input, output);

   { Test n-dimensional Runge Kutta Fehlberg integrator. }

{$include:'rkf45.pas'}

var
   x_init,         { initial value of independent variable }
   x_fin,          { final   value of independent variable }
   x_low,          { low  intermediate value of indep. variable }
   x_high          { high intermediate value of indep. variable }
      : real;
   y_init,         { initial value of dep. variable }
   y_low,          { initial value of dep. variable for step }
   y_high,         { value of dep. var. after step }
   y_fin,          { final answer, if integration was successful }
   y_prime,        { storage space for derivatives }
   abs_err,        { absolute error tolerated on integration }
   rel_err,        { relative error tolerated on integration }
   epsilon         { total tolerance on integration }
      : workarray;
   step,           { step size given to rkf45 }
   new_step,       { step size recommended by rkf45 }
   min_step,       { smallest tolerable step size }
   k               { parameter for f() }
      : real;  
   finished        { status of stepping }
      : boolean;
   n,              { number of dimensions }
   fail,           { status of rkf45 }
   m,              { used for computing starting step size }
   i               { dummy index for array indexing }
      : integer;

function f2
   (n : integer; x : real;   y : workarray; var y_prime : workarray)
      :  boolean;
  { Generate exponential function from first order eqn.

     y' := k  y

  }

{shared
   k : real;           constant for function
}

begin
   if (k * x < 230.0)  then { within good range of x }
   begin
      y_prime[1] := k * y[1];
      f2 := true;
   end
   else
      f2 := false;  { say solution blows up }
end { end f2 };

function f3
   (n : integer; x : real;   y : workarray; var y_prime : workarray)
      :  boolean;
  { Generate trig. function from second order eqn.

              2
     y'' := -k   y

   by transforming to 2 first order equations,

     y'  := y
      1      2

              2
     y'  := -k  y
      2          1
  }

{shared
   k : real;           constant for function
}

begin
   if (k * x < 230.0)  then { within good range of x for demo }
   begin
      y_prime[1] := y[2];
      y_prime[2] := -k * k * y[1];
      f3 := true;
   end
   else
      f3 := false;  { say solution blows up  for demo. }
end { end f3 };

begin      { main }
   { Initial conditions }
   { n := 1;             single first  order equation for f2 }
   n := 2;             { single second order equation for f3 }

   k := 6.283185308;
   x_init := 0.0;
   y_init[1] := 1.0;   { y(0) }
   y_init[2] := 0.0;   { generate cosine(2 pi x) }
   x_fin  := 10.0;     { integrate to x = 10 }

   { Initialization for integration itself }
   x_high := x_init;
   for i := 1 to n do
      y_high[i] := y_init[i];
   m :=  25;
   new_step := (x_fin - x_init) / m;
   rel_err[1] := 1.0e-3;
   abs_err[1] := 1.0e-3;
   rel_err[2] := k * rel_err[1];  { Must allow different tolerance }
   abs_err[2] := k * abs_err[1];  { on the derivative. }
   min_step := 1.0e-5;
   fail := 0;
   finished := false;

   { Stepping loop with automatic step size selection. }
   while ((fail = 0)  and not  finished) do
   begin
      step := new_step;
      x_low := x_high;
      x_high := x_low + step;
      for i := 1 to n do
      begin
         y_low[i] := y_high[i];
         epsilon[i] := abs_err[i] + abs(y_low[i]) * rel_err[i];
      end;

      if (((x_high - x_fin) / step) > 0) then
      begin
         x_high := x_fin;
         finished := true;
      end;
      rkf45(n, f3, y_prime, x_low, y_low, x_high, epsilon, step,
             min_step, y_high, new_step, fail);

      { to record progress of integration  activate the following: 
       
         writeln(' ', (x_high - x_low):14, ' ', x_high:8:3, ' ',
                 y_high[1]:14, ' ', y_high[2]:14, ' ', new_step:14);   }
   end; { end stepping } ;

   if (fail = 0) then   { integration was successful }
   begin
      writeln(' Solution = ');
      for i := 1 to n do
      begin
         y_fin[i] := y_high[i];
         write(' ', y_fin[i]:14, ' ');
      end;
      writeln;
   end
   else
      writeln(' RKF45 failed to find a solution');

end  { end main } .