function [y, t] = eulerSys (f, t0, tlast, y0, N)
% function [y, t] = eulerSys (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 1st-order Euler 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'
%             It is necessary f accept inputs
%               t  scalar (time)
%               y  vector (state)
%             and output a vector of the same size and shape as 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       matrix of approximate value of solution at mesh points
%             The jth column corresponds to the jth time.
%     t       mesh points

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