function retval = ClassERK(fName, t0, y0, hStep, nSteps, fEvery)
% usage: retval = ClassERK(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 classical explicit
% Runge-Kutta third order.  The optional flag fEvery is default
% set to false.  User can change it via parameter list. Return value
% in this case is an n x (nsteps+1) matrix whose columns are 
% computed solutions interpolated to the grid of times 
% t = t0 + k*hStep, k = 0,...,nSteps+1. 

% programmer: t. shores
% last rev: 4/10/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
% t: current time value
% y: most recent value of computed solution
% jj: index variable
% ff: computed intermediate function values in ERK algorithm
% fEvery: flag to return initial and every subsequent profile in a matrix

global gError; 

% initialization code
y0 = y0(:);   % ensure column
n = size(y0, 1);
if (exist('fEvery', 'var'))
  retval = zeros(n, nSteps + 1);
  retval(:,1) = y0(:);
  fEvery = true;
else
  fEvery = false;
end
c = [0; 0.5; 1]; 
b = [1/6; 2/3; 1/6];
a = [0,0,0;0.5,0,0;-1,2,0];
ff = zeros(n,3);
t = t0;
y = y0;
h = hStep;
% main loopy
for jj = 1:nSteps
  % do the ERK3 step
  ff(:,1) = feval(fName, t, y);
  ff(:, 2) = feval(fName, t + c(2)*hStep, y + hStep*a(2,1)*ff(:, 1));
  ff(:, 3) = feval(fName, t + c(3)*hStep, y + hStep*(a(3,1)*ff(:, 1) + a(3,2)*ff(:, 2)));
  y = y + hStep*(b(1)*ff(:, 1) + b(2)*ff(:, 2) + b(3)*ff(:, 3));
  t = t + hStep;
  if (fEvery)
    retval(:, jj + 1) = y;
  end
end;
if (~fEvery)
  retval = y;
end

end
