function [y, t] = heun(f, t0, tlast, y0, N)
% function [y, t] = heun(f, t0, tlast, y0, N)
%
%  This routine calculates an approximate solution at mesh points
%  t0, t1, ..., tN  to the initial-value problem
%
%          y' = f(t, y),   subject to y(t0) = y0,
%
%  using the 2nd-order Runge-Kutta method also known as Heun's
%  method.  We assume evenly-spaced mesh points, with stepsize
%
%              h = (tlast - t0) / N.
%
%   INPUTS:
%     f       function handle to the RHS in the ODE giving y'
%     t0      initial "time"
%     tlast   last "time" at which to approximate the solution
%     y0      initial value, the value of the solution at the initial time
%     N       number of steps to take to go from t0 to tN
%
%  OUTPUTS:
%     y       approximate value of solution at mesh points
%     t       mesh points

  h = (tlast - t0) / N;
  t = [t0:h:tlast]';
  y = y0;
  for jj = 1:N
    k1 = f(t(jj), y(jj));
    k2 = f(t(jj + 1), y(jj) + h*k1);
    y(jj + 1) = y(jj) + h*(k1 + k2) / 2;
  end
  y = y(:);
end