function retval = ERK23(fName, t0, y0, hStep, nSteps, fEvery)
% usage: retval = ERK23(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 1 vector and the method of solution is explicit embedded
% third order Runge-Kutta.  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/10/16

% Octave and Matlab compatible file
% user defined local variables:
% s: factor for downsizing variable stepsize (default s = 1/2)
% S: factor for upsizing variable stepsize (default S = 2)
% eMax: maximum error per unit step tolerance (default eMax=0.01)
% eMin: minimum error per unit stepsize tolerance (default eMin=0.1*eMax)
% hMin: minimum acceptable stepsize (default hMin = .0001)
% remaining local variables:
% gError: signal for failure of routine 
% n: rows of y0
% nextNode:  next node to be interpolated if fEvery is true
% nextIndex: index of next node to be interpolated
% aa: RK matrix -- double letters mean full embedded tableau
% b, bb: RK weights
% cc: RK nodes
% ff: computed intermediate function values in ERK algorithm
% t: current time value to which y values are accepted
% y: current value of y which is accepted at time t
% h: current stepsize 
% tLast: finishing time for this function
% ii: index variable
% d1, d2, d3: divided differences for interpolation

global gError; 

% user defined local variables
s = 0.5;
S = 2;
eMax = 1e-4;
eMin = 0.1*eMax;
hMin = 0.001*eMax;
% remaining local variables
myEps = 10*eps;
gError = false;
n = size(y0, 1);
if (exist('fEvery', 'var'))
  retval = zeros(n, nSteps + 1);
  retval(:,1) = y0(:);
  nextNode = t0 + hStep;
  nextIndex = 2;
  fEvery = true;
else
  fEvery = false;
end
aa = [0, 0, 0; 2/3, 0, 0; 0, 2/3, 0]; % really atilde. 
% this is a very simple ERK23 -- Bogati-Shampine is Matlab's ERK23 method
b = [1/4; 3/4]; % really b
bb = [1/4; 3/8; 3/8];  % really btilde, and fyi since not needed
cc = [0;2/3;2/3]; % really ctilde
ff = zeros(n,3);
t = t0;
tLast = t0 + hStep*nSteps;
y = y0(:);  % make sure y is a column
h = hStep;
nextNode = t0 + hStep;
% main loop
while (t < tLast)
  % calculate next y
  ff(:,1) = feval(fName, t, y);
  ff(:,2) = feval(fName, t + cc(2)*h, y + aa(2, 1)*h*ff(:,1));
  ff(:,3) = feval(fName, t + cc(3)*h, y + aa(3, 2)*h*ff(:,2));
  yNew = y + h*(b(1)*ff(:, 1) + b(2)*ff(:, 2));
  tNew = t + h;
  % calculate local error (this is wired to current tableau values)
  kappa = 3*h*norm(ff(:,2)-ff(:,3),inf)/8;
  if (kappa < h*eMax)
    % error is acceptable, so prepare to place values in Nodes arrays
    if (fEvery)
      % but first we need to do all necessary back interpolation
      while((nextNode <= tNew)  && (nextIndex <= nSteps + 1))
        % back up for interpolation to nextNode
        hBak = h;
        h = tNew - nextNode;
        % calculate next retval
        ff(:,1) = feval(fName, t, y);
        ff(:,2) = feval(fName, t + cc(2)*h, y + aa(2, 1)*h*ff(:,1));
        retval(:, nextIndex) = y + h*(b(1)*ff(:, 1) + b(2)*ff(:, 2));
        nextNode = nextNode + hStep;
        nextIndex = nextIndex + 1;
        h = hBak;
      end   % <= tNew
      if (nextIndex > nSteps + 1)
        break;  % done, so check out
      end;
    end % fEvery
    % check for finished condition
    if (tLast <= tNew)
      % done, so interpolate and check out
       h = tNew - tLast;
       ff(:,1) = feval(fName, t, y);
       ff(:,2) = feval(fName, t + cc(2)*h, y + aa(2, 1)*h*ff(:,1));
       y = y + h*(b(1)*ff(:, 1) + b(2)*ff(:, 2));
       break;
    end
    % if we're here, accept the updates, but we're not done
    t = t + h;
    y = yNew;
    if (kappa < h*eMin)
      % error is too small so size up h for next round
      h = S*h;
    end % kappa < h*eMin  
    if (h < hMin)
      % stepsize too small, so resize and break out of the while
      gError = true;
    end; % h < hMin
    % if we get here, new t and h are accepted
  else % end acceptable error
    h = s*h;  % error is unacceptable
  end; 
end;  % end while
% output section
% if we ran into an error condition, print a message
if (gError)
  disp('Error 2 in RK23: stepsize too small at time');
  disp(t)
  retval = y;
  return;
end;
if (~fEvery)
  retval = y;
end

end
