function [t,y] = rk4ode(dydt,tspan,y0,h,varargin)
%
% rk4ode: 4th order Runge-Kutta ODE solver 
%   [t,y] = rk4ode(dydt,tspan,y0,h,p1,p2,...):
%           uses 4th order RK method to integrate an ODE
%
% input:
%   dydt = function handle that evaluates the ODE 
%   tspan = [ti, tf] where ti and tf = initial and
%           final values of independent variable
%   y0 = initial value of dependent variable
%   h = step size
%   p1,p2,... = additional parameters used by dydt
%
% output:
%   t = vector of independent variable
%   y = vector of solution for dependent variable
%

if nargin<4,error('at least 4 input arguments required'),end

ti = tspan(1);tf = tspan(2);
if ~(tf>ti),error('upper limit must be greater than lower'),end

t = ti:h:tf; 
n = length(t);

% if necessary, add an additional value of t
% so that range goes from t = ti to tf
if t(n)<tf  
  t(n+1) = tf;
  n = n+1;
end

y = y0; % initialize y
for i = 1:n-1  % implement RK4 method
  k1 = dydt(t(i),y(i),varargin{:});  
  ymid = y(i)+k1*h/2;
  k2 = dydt(t(i)+h/2,ymid,varargin{:}); 
  ymid = y(i)+k2*h/2;
  k3 = dydt(t(i)+h/2,ymid,varargin{:}); 
  yend = y(i)+k3*h;
  k4 = dydt(t(i)+h,yend,varargin{:}); 
  phi = (k1+2*k2+2*k3+k4)/6;
  y(i+1) = y(i) + phi*(t(i+1)-t(i));
end

