function [ x, y ] = rk3 ( f_ode, x_range, y_initial, numSteps )
% usage:  [ x, y ] = rk3 ( f_ode, x_range, y_initial, numSteps )
% Solve a system of ODE's using a third-order Runge-Kutta method.
% f_ode = function handle to compute the right side of the ODE
% x_range   = vector of length 2 describing an interval over
%             which the ODE is to be solved
% y_initial = column vector of initial values
% numSteps     = number of intervals
% x         = vector of length (numSteps+1) spanning the range given
%             by t_range
% y         = matrix with one row for each value of x, the 
%             contents of the row are the solution at time x.

% M. M. Sussman
% $Date: 2014/04/09 14:29:58 $

% make sure that x_range is precisely two values
if length(x_range) ~= 2
  x_range
  error('RK3: x_range must be a vector with exactly two entries.')
elseif x_range(1)>=x_range(2)
  x_range
  error('RK3: x_range(2) must be larger than x_range(1).')
end

% make sure numSteps is positive integer
if (round(numSteps) ~= numSteps) | (numSteps <= 0)
  numSteps
  error('RK3: numSteps must be a positive integer.')
end

% check that y_initial is a column vector
s=size(y_initial);
if s(1) > s(2)
  error('RK3: y_initial must be a column vector')
end

% check that f_ode returns a column vector of the correct size
rhs=f_ode(x_range(1),y_initial);
if size(rhs) ~= size(y_initial)
  error('RK3: right side function must return a column vector.')
end

x=zeros(1,numSteps+1);
y=zeros(size(y_initial,1),numSteps+1);
x(1) = x_range(1);
h = ( x_range(2) - x_range(1) ) / numSteps;
y(:,1) = y_initial;
for k = 1 : numSteps
  xa = x(1,k)+h/2 ;
  ya = y(:,k)+.5*h*f_ode(x(1,k),y(:,k)) ;

  xb = x(1,k)+h;
  yb = y(:,k)+h*(2*f_ode(xa,ya)-f_ode(x(1,k),y(:,k)));

  x(1,k+1) = x(1,k) + h;
  y(:,k+1) = y(:,k) + h/6 * (f_ode(x(1,k),y(:,k))+4*f_ode(xa,ya)+f_ode(xb,yb));
end
