function retval = ClassERK(fName, t0, y0, hStep, nSteps)
% usage: retval = ClassERK(fName,t0,y0,hStep,nSteps)
% 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 internal flag fEvery is default set
% to false.  User can change it, in which case the returned value is
% an mX(nsteps+1) matrix whose columns are y0 followed by the
% solution at every computed step.

% programmer: t. shores
% last rev: 4/4/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: lower triangle of RK matrix of k coefficients in RK stages
% t: current time value
% jj: index variable
% kCoefs: function evaluations for the method
% fEvery: flag to return initial and every subsequent profile in a matrix
% newVal: most recent value of computed solution

global gError; 

% initialization code
fEvery = false;  % set to true for full profile
c = [0; 0.5; 1]; 
b = [1/6; 2/3; 1/6];
a = [0,0,0;0.5,0,0;-1,2,0];
kCoef = zeros(length(y0), 3);
t = t0;
newVal = y0(:); % make y0 a column
if (fEvery)
  retval = zeros(length(y0), nSteps + 1);
end
% main loop
for jj = 1:nSteps
  % do the ERK3 step
  kCoef(:,1) = feval(fName, t, newVal);
  kCoef(:, 2) = feval(fName, t + c(2)*hStep, newVal + hStep*a(2,1)*kCoef(:, 1));
  kCoef(:, 3) = feval(fName, t + c(3)*hStep, newVal + hStep*(a(3,1)*kCoef(:, 1) + a(3,2)*kCoef(:, 2)));
  newVal = newVal + hStep*(b(1)*kCoef(:, 1) + b(2)*kCoef(:, 2) + b(3)*kCoef(:, 3));
  t = t + hStep;
  if (fEvery)
    retval(:, jj + 1) = newVal;
  end
end;
if (~fEvery)
  retval = newVal;
end

end
