function [t,y]=Heun(f,tspan,y0,n);
% HEUN solves system of ODEs using Heun's method
%   [t,y]=Heun(f,tspan,y0,n) solves dy/dt=f(t,y) with initial
%   value y0 on the time interval tspan doing n steps of Heun's 
%   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)+dt/3,y(:,i)+dt/3*k1);
  k3=f(t(i)+2*dt/3,y(:,i)+2*dt/3*k2);
  y(:,i+1)=y(:,i)+dt*(1/4*k1+3/4*k3);
end;
t=t';y=y';                        % to return results in columns

