function [y, t] = rk4sys(f, tspan, y0, N)
% function [y, t] = rk4sys(f, tspan, y0, N)
%
%  This routine calculates an approximate solution at mesh points
%  t0, t1, ..., tN  to the vector initial-value problem
%
%          y' = f(t, y),   subject to y(t0) = y0,
%
%  using the most familiar 4th-order Runge-Kutta method.
%  We assume evenly-spaced mesh points, with stepsize
%
%              h = (tspan(2) - tspan(1)) / 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.
%     tspan   vector containing initial and final "time"
%     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

  t0 = tspan(1); tlast = tspan(2);
  h = (tlast - t0) / N;
  t = [t0:h:tlast]';
  y = y0(:);       % ensure initial vector is row vector
  for jj = 1:N
    k1 = f(t(jj), y(:, jj));
    k2 = f(t(jj) + h/2, y(:, jj) + h*k1/2);
    k3 = f(t(jj) + h/2, y(:, jj) + h*k2/2);
    k4 = f(t(jj + 1), y(:, jj) + h*k3);
    y(:, jj + 1) = y(:, jj) + h*(k1 + 2*k2 + 2*k3 + k4) / 6;
  end
end