-- rk4th.ada  RUNGE-KUTTA 4th order integration of differential equation
--            given a function, dy/dx = y' = f(x,y)
--            given an initial point x0, y0 and step size h and number steps s
--            determine the final point xf, yf
generic
  type F_to_int is access function (X , Y : Float) return Float; 
procedure RK4thgp(F : F_to_int;     -- function to integrate
                  X : in out float; -- initial and final value
                  Y : in out float; -- initial and final value
                  H : float;        -- step size
                  S : integer);     -- number of steps

procedure RK4thgp(F : F_to_int;     -- function to integrate
                  X : in out float; -- initial and final value
                  Y : in out float; -- initial and final value
                  H : float;        -- step size
                  S : integer) is   -- number of steps
  K1, K2, K3, K4 : Float;
begin
  for I in 1..S loop
    K1 := H*F(X,Y);
    K2 := H*F(X+H/2.0,Y+K1/2.0);
    K3 := H*F(X+H/2.0,Y+K2/2.0);
    K4 := H*F(X+H,Y+K3);
    Y  := Y+(K1+2.0*K2+2.0*K3+K4)/6.0;
    X  := X+H;
  end loop;
end RK4thgp;


with RK4thgp;              -- with the generic procedure
with Text_IO; use Text_IO;
procedure RK4thg is
  X, Y, H : Float;
  S : Integer;
  type F_2_f is access function (X , Y : Float) return Float; 
  procedure RK4_Integrate is new RK4thgp(F_2_f);


  function T1(X, Y : Float) return Float is
  begin
    return 2.0*Y/X;   -- initial condition x=1, y=2
                      -- solution y = 2x**2
  end T1;


  function T2(X, Y : Float) return Float is
  begin
    return Y*(Y*Y+3.0*X*X)/(2.0*(X**3)); -- initial condition x=1/2, y=1/2
    -- solution y**2 = (x**3)/(1-x) 
  end T2;
begin
  Put_Line("Runge-Kutta solution");
  X := 1.0;
  Y := 2.0;
  H := 0.1;
  S := 10;
  RK4_Integrate(T1'access, X, Y, H, S);
  Put_Line("Final X = " & Float'Image(X) & ", Y = " & Float'Image(Y));
  Put_Line("Exact X = 2.0, Y=8.0"); New_Line;
  X := 0.5;
  Y := 0.5;
  H := -0.025; -- X gets smaller
  S := 10;
  RK4_Integrate(T2'access, X, Y, H, S);
  Put_Line("Final X = " & Float'Image(X) & ", Y = " & Float'Image(Y));
  Put_Line("Exact X = 0.25, Y=0.1443375    1/(4*sqrt(3))");
end RK4thg;