function [t,y]=HeunNegative(f,tspan,y0,n);
% HEUNNEGATIVE solves system of ODEs using Heun's 2nd possible 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 
%   2nd possible method with negative weights

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=feval(f,t(i),y(:,i));
  k2=feval(f,t(i)+dt/3,y(:,i)+dt/3*k1);
  k3=feval(f,t(i)-dt/3,y(:,i)-dt*k1+2*dt/3*k2);
  y(:,i+1)=y(:,i)+dt*(-2*k1+9/4*k2+3/4*k3);
end;
t=t';y=y';                        % to return results in columns

