function [tt,u]=rk4(odefun,tspan,y0,Nh,varargin);
%RK4 Solve differential equations using the Runge-Kutta
% 4th order explicit scheme
%  [T,Y]=RK4(ODEFUN,TSPAN,Y0,NH) with TSPAN=[T0,TF]
%  integrates the system of differential equations
%  y'=f(t,y) from time T0 to TF with initial condition
%  Y0 using RK4 method on an equispaced
%  grid of NH intervals.Function ODEFUN(T,Y) must return
%  a column vector corresponding to f(t,y). Each row in
%  the solution array Y corresponds to a time returned
%  in the  column vector T.
%  [T,Y] = RK4(ODEFUN,TSPAN,Y0,NH,P1,P2,...) passes
%  the additional parameters P1,P2,... to the function
%  ODEFUN as ODEFUN(T,Y,P1,P2...).

tt=linspace(tspan(1),tspan(2),Nh+1);
h=(tspan(2)-tspan(1))/Nh; hh=h*0.5; 
y=y0(:);
w=y; u=y.';
h6=h/6;
for t=tt(1:end-1)
  
  k1=odefun(t,w,varargin{:});
  t1 = t + hh; w1 = w + hh* k1;
  k2=odefun(t1,w1,varargin{:});
  w1 = w + hh* k2;
  k3=odefun(t1,w1,varargin{:});
  t1 = t + h; w1 = w + h*k3;
  k4=odefun(t1,w1,varargin{:});
  w=w + h6*(k1+2*k2+2*k3+k4);
  u = [u; w.'];

end

tt=tt';

