function [T, X] = rkf( f, a, b, x0, tol, hmax, hmin )
  % --------------------------------------------------------------------------
  % Jonathan R. Senning <jonathan.senning@gordon.edu>
  % Gordon College
  % March 23, 1999
  %
  % $Id: rkf.m,v 1.1 2005/04/27 12:58:44 senning Exp senning $
  %
  % Usage: [t, x] = rkf(@f, a, b, x0, tol, hmax, hmin)
  %
  % Returns:
  %	t:	array of independent variable values (irregular steps sizes)
  %	x:	array of corresponding solution function values
  %
  % Parameters:
  %	f:	handle of function enclosed in quotes equal to dx/dt = f(x,t)
  %	a:	left-hand endpoint of interval (initial condition is here)
  %	b:	right-hand endpoint of interval
  %	x0:	initial x value: x0 = x(a)
  %	tol:	maximum value of local truncation error estimate
  %	hmax:	maximum step size
  %	hmin:	minimum step size
  %
  % This function implements 4th-5th order Runge-Kutta-Fehlberg Method to
  % solve the initial value problem
  %
  %       dx
  %       -- = f(x,t),     x(a) = x0
  %       dt
  % 
  % on the interval [a,b].
  %
  % Based on pseudocode presented in "Numerical Analysis", 6th Edition, by
  % Burden and Faires, Brooks-Cole, 1997.
  % --------------------------------------------------------------------------

  % Coefficients used to compute first argument of f

  a2  =    1 /  4;
  a3  =    3 /  8;
  a4  =   12 / 13;
  a5  =    1;
  a6  =    1 /  2;

  % Coefficients used to compute second argument of f

  b21 =      1 /    4;
  b31 =      3 /   32;
  b32 =      9 /   32;
  b41 =   1932 / 2197;
  b42 =  -7200 / 2197;
  b43 =   7296 / 2197;
  b51 =    439 /  216;
  b52 =     -8;
  b53 =   3680 /  513;
  b54 =   -845 / 4104;
  b61 =     -8 / 27;
  b62 =      2;
  b63 =  -3544 / 2565;
  b64 =   1859 / 4104;
  b65 =    -11 /   40;

  % Coefficients used to compute local truncation error estimate.  These
  % come from subtracting a 4th order RK estimate from a 5th order RK estimate.

  r1  =      1 /  360;
  r3  =   -128 /  4275;
  r4  =  -2197 / 75240;
  r5  =      1 /    50;
  r6  =      2 /    55;

  % Coefficients used to compute 4th order RK estimate

  x1  =     25 /  216;
  x3  =   1408 / 2565;
  x4  =   2197 / 4104;
  x5  =     -1 /    5;

  % Define range for stepsize factor

  deltamin = 0.1;
  deltamax = 4.0;
  
  % Define boolean values

  false = 0;
  true = ~false;	% means: True is equal to "NOT False"

  % Set t and x according to initial condition and assume that h starts
  % with a value that is as large as possible.

  t = a;
  x = x0;
  h = hmax;
  done = false;

  % Initialize arrays that will be returned

  T = t;
  X = x;

  while ( ~done )
%    printf( "t = %f, h = %e\n", t, h ); fflush(stdout);

    % Compute values needed to compute truncation error estimate and
    % the 4th order RK estimate.

    k1 = h * f( x, t );

    t0 = t + a2 * h;
    x0 = x + b21 * k1;
    k2 = h * f( x0, t0 );

    t0 = t + a3 * h;
    x0 = x + b31 * k1 + b32 * k2;
    k3 = h * f( x0, t0 );

    t0 = t + a4 * h;
    x0 = x + b41 * k1 + b42 * k2 + b43 * k3;
    k4 = h * f( x0, t0 );

    t0 = t + a5 * h;
    x0 = x + b51 * k1 + b52 * k2 + b53 * k3 + b54 * k4;
    k5 = h * f( x0, t0 );

    t0 = t + a6 * h;
    x0 = x + b61 * k1 + b62 * k2 + b63 * k3 + b64 * k4 + b65 * k5;
    k6 = h * f( x0, t0 );

    % Compute the estimate of the local truncation error.  If it's small
    % enough then we accept this step and save the 4th order estimate.

    printf("%e %e %e %e %e\n", k1, k3, k4, k5, k6 );
    r = abs( r1 * k1 + r3 * k3 + r4 * k4 + r5 * k5 + r6 * k6 ) / h;
    if ( r <= tol )
      t = t + h;
      x = x + x1 * k1 + x3 * k3 + x4 * k4 + x5 * k5;
      T = [T, t];
      X = [X, x];
    end

    % Now compute next step size, and make sure that it is not too big.

%    printf("tol = %e, r = %e\n", tol, r );
    delta = ( 0.5 * tol / r )^0.25;

    if ( delta <= deltamin )
      delta = deltamin;
    elseif ( delta >= deltamax )
      delta = deltamax;
    end

    h = delta * h;

    if ( h > hmax )
      h = hmax;
    end

%    printf( "delta = %e, h = %e\n", delta, h );

    % Check to see if we are done, or within a single step of being done.
    % If not, then make sure that h is not too small -- if it is then we
    % terminate with an error message.

    if ( t >= b )
      done = true;
    elseif ( t + h > b )
      h = b - t;
    elseif ( h < hmin )
      done = true;
      fprintf( 'h is smaller than %e: function terminated.\n', hmin );
    end

  end

end
