function retval = ERK4(fName, t0, y0, hStep, nSteps, fEvery)
% usage: retval = ERK4(fName, t0, y0, hStep, nSteps, fEvery))
% description: Returns an approximate solution to dy/dt = fcn(t,y) at 
% t = t0 + hStep*nSteps.  Here fName is the name of fcn, y(t0) = y0
% an n x 1column vector and the method of solution is explicit
% Runge-Kutta fourth order.  The optional flag fEvery is default
% set to false.  User can change it via parameter list. Return value
% in this case is an nX(nsteps+1) matrix whose columns are 
% computed solutions interpolated to the grid of times 
% t = t0 + k*hStep, k = 0,...,nSteps+1. In case of an error condition, 
% an error message is displayed and the global gError is set to true. 

% programmer: t. shores
% last rev:4/7/16

% Octave and Matlab compatible file
% local variables
% a, b, c, comprise the Butcher tableau for this method
% c: column vector of time nodes (time coefficients in function evaluation)
% b: row vector of weights (function coefficients in RK sum)
% a: RK matrix of k coefficients in RK stages
% jj: index variable
% ff: function evaluations for the method
% gError: flag to return initial and every subsequent profile in a matrix
% t: current time value to which y values are accepted
% y: current value of y which is accepted at time t
% h: current stepsize 

global gError; 

% initialization code
n = size(y0, 1);
if (exist('fEvery', 'var'))
  retval = zeros(n, nSteps + 1);
  retval(:,1) = y0(:);
  fEvery = true;
else
  fEvery = false;
end
% coefficients for 3/8 variation on RK4
c = [1/2; 1/2; 1];      % c = [1/3; 2/3; 1];          % c = [1/2; 1/2; 1];      % for plain RK4
b = [1/6, 1/3, 1/3, 1/6]; % b =  [1/8, 3/8, 3/8, 1/8];   % b = [1/6, 1/3, 1/3, 1/6]; % for plain RK4
a = [1/2, 0, 0; 0, 1/2, 0; 0, 0, 1]; % a = [1/3, 0, 0; -1/3, 1, 0; 1, -1, 1];
%a = [1/2, 0, 0; 0, 1/2, 0; 0, 0, 1]; % lower triangle of plain RK4
ff = zeros(length(y0),4);
t = t0;
y = y0(:);  % make sure y is a column
h = hStep;
% main loop
for jj = 1:nSteps
  % do the ERK4 step
  ff(:,1) = feval(fName, t, y);
  ff(:, 2) = feval(fName, t + c(1)*hStep, y + hStep*a(1,1)*ff(:, 1));
  ff(:, 3) = feval(fName, t + c(2)*hStep, y + hStep*(a(2,1)*ff(:, 1) + a(2,2)*ff(:, 2)));
  ff(:,4) = feval(fName, t + c(3)*hStep, y + hStep*(a(3,1)*ff(:, 1) + (a(3,2)*ff(:, 2) + a(3,3)*ff(:, 3))));
  y = y + hStep*(b(1)*ff(:, 1) + b(2)*ff(:, 2) + b(3)*ff(:, 3) + b(4)*ff(:, 4));
  t = t + hStep;
  %  do fill-in if requested
  if (fEvery)
     retval(:, jj + 1) =y;
  end
end
if (~fEvery)
  retval = y;
end

end
