function x = rk4( f, x0, t )
  % --------------------------------------------------------------------------
  % Jonathan R. Senning <jonathan.senning@gordon.edu>
  % Gordon College
  % March 31, 1999
  %
  % Usage: x = rk4(@f, x0, t)
  %
  % Returns:
  %     x:      array of solution function values corresponding to the
  %             values in the supplied array t.  Note that x will have
  %             the same row dimension as x0 (number of equations) and
  %             the same column dimesion as t (number of points).
  %
  % Parameters:
  %     f:      handle of function dx/dt = f(x,t).  This function should
  %             return a quantity that has the same dimensions as x,
  %             namely Mx1 where M is the number of equations.
  %     x0:     initial x value: x0 = x(t(1)).  Notice that if a system of
  %             equations is being solved then this should be a column vector
  %             (dimensioned Mx1 where M is the number of equations).
  %     t:      array of time values to determine x at.  This array should
  %             be dimesioned 1xN where N is the number of values in t.
  %
  % This function implements 4th order Runge-Kutta to solve
  %
  %       dx
  %       -- = f(x,t),     x(t(1)) = x0
  %       dt
  %
  % at the t values stored in the t array (so the interval of solution is
  % [t(1), t(N)]. Notice that it works equally well for scalar functions
  % f(x,t) (in the case of a single 1st order ODE) or for vector functions
  % f(x,t) (in the case of multiple 1st order ODEs).
  % --------------------------------------------------------------------------

  neqn = size( x0, 1 );
  n = length( t );
  x = zeros( neqn, n );
  x(:,1) = x0;

  for i = 1 : n-1
    h = t(i+1) - t(i);
    k1 = h * f( x(:,i), t(i) );
    k2 = h * f( x(:,i) + 0.5 * k1, t(i) + 0.5 * h );
    k3 = h * f( x(:,i) + 0.5 * k2, t(i) + 0.5 * h );
    k4 = h * f( x(:,i) + k3, t(i+1) );
    x(:,i+1) = x(:,i) + ( k1 + 2.0 * ( k2 + k3 ) + k4 ) / 6.0;
  end

end
