function [t,x] = heun(xprime, tspan, x0, h)

% function [t,x] = heun(xprime, [t0 tfinal], x0, h)
% 
% Approximate the solution to the ODE:  x'(t) = xprime(x,t)
% from t=t0 to t=tfinal with initial condition x(t0)=x0
% using Heun's method (improved Euler's method).
% xprime should be a function that can be called like: xprime(t,x).
% h is the stepsize:  reduce h to improve the accuracy.
% The ODE can be a scalar or vector equation.

 t0 = tspan(1); tfinal = tspan(end);

% set up the t values at which we will approximate the solution
 t=[t0:h:tfinal];

% include tfinal even if h does not evenly divide tfinal-t0
 if t(end)~=tfinal, t = [t tfinal]; end    

% execute Heun's method
 x = [x0 zeros(length(x0),length(t)-1)];
 for j=1:length(t)-1
    k1 = feval(xprime,t(j),x(:,j));
    x(:,j+1) = x(:,j) + .5*h*(k1 + feval(xprime,t(j)+h,x(:,j)+h*k1)); 
 end