function [t,y]=RK4(f,tspan,y0,n);
% RK4 solves system of ODEs using classical 4th order Runge Kutta
%   [t,y]=RK4(f,tspan,y0,dt) solves dy/dt=f(t,y) with initial
%   value y0 on the time interval tspan doing n steps of the 
%   classical 4th order Runge Kutta method.

dt=(tspan(2)-tspan(1))/n; 
t=tspan(1):dt:tspan(2);   
y(:,1)=y0(:);                       % colon to make column vector
for i=1:n,
  k1=f(t(i),y(:,i));
  k2=f(t(i)+0.5*dt,y(:,i)+0.5*dt*k1);
  k3=f(t(i)+0.5*dt,y(:,i)+0.5*dt*k2);
  k4=f(t(i)+dt,y(:,i)+dt*k3);
  y(:,i+1)=y(:,i)+dt*((k1+k4)/6+(k2+k3)/3);
end;
t=t';y=y';                          % to return results in columns