% function y = heun(a,b,n,y0,f)
% 
%  uses Heun's method to compute solution to IVP
%
%    y'    = f(t,y),  a <= t <= b
%    y(a)  = y0
%
% inputs:
%  a,b   time interval to compute solution
%  n     number of time steps to take inside [a,b]
%  y0    initial condition at t=a
%  f     function handle, for example
%
%        f = @(t,y) t^2 + y/2;
% outputs:
%  y     vector of size n+1 y(i) contains the approximation to
%        y(a + (i-1)*h), for i=1,..,n+1, with h=(b-a)/n
function y = heun(a,b,n,y0,f)
 y(1) = y0;
 h=(b-a)/n;
 for i=1:n,
  ti=a+(i-1)*h;
  % TO BE COMPLETED
  y(i+1) = y(i) +  % TO BE COMPLETED
 end;
