function [t,y] = heun(f,tspan,y0,h)
%HEUN solve differential equation with Heun's method 
%   [T, Y] = HEUN(F,TSPAN,Y0,H) with TSPAN = [T0, TFINAL] integrates
%   the system of differential equations y' = f(t,y) from time T0 to
%   TFINAL with initial conditions Y0 employing a constant time step
%   H. F is a function handle. For a scalar T and a vector Y, F(T,Y)
%   must return a column vector corresponding to f(t,y). Each row in
%   the solution array Y corresponds to a time in the vector T.
