function y = HeunODE(f,x,y0)
%
% HeunODE uses Heun's method to solve a first-order ODE given in the form
% y' = f(x,y) subject to initial condition y0.
%
%  y = HeunODE(f,x,y0) where
%
%       f is an inline function representing f(x,y),
%       x is a vector representing the mesh points,
%       y0 is a scalar representing the initial value of y,  
%
%       y is the vector of solution estimates at the mesh points.

y = 0*x;    % Pre-allocate
y(1) = y0; h = x(2)-x(1);
for n = 1:length(x)-1,
    k1 = f(x(n),y(n));
    k2 = f(x(n)+h,y(n)+h*k1);
   y(n+1) = y(n)+h*(k1+k2)/2;
end