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/7/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
% tNodes: array of two most recent accepted time values
% yNodes: array of  two most recent accepted time values
% 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
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;
tNodes = zeros(1, 2);
tNodes(1,2) =  t0;
tNodes(1,1) = t0 - hStep; 	% artificial time for possible quadratic interpolation
tLast = t0 + hStep*nSteps;
y = y0(:);  % make sure y is a column
yNodes = zeros(n, 2);
yNodes(:, 2) = y;
yNodes(:, 1) = y;	% artificial node 
h = 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));
  % 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
    t = t + h;
    y = yNew;
    % but first we need to do all necessary back interpolation
    if (fEvery)
       while (nextNode <= t)
          % ok, interpolation time  
          d1 = (yNodes(:, 2) - yNodes(:, 1))/(tNodes(2) - tNodes(1));
          d2 = (y - yNodes(:, 2))/(t - tNodes(2));
          d3 = (d2 - d1)/(t - tNodes(1));
          retval(:, nextIndex) = yNodes(:, 1) + d1*(nextNode - tNodes(1) )+ d3*(nextNode - tNodes(1))*(nextNode - tNodes(2));
          % check for finished condition to break out of the while
          if (t >= tLast)
            break;
          end
          nextNode = nextNode + hStep;
          nextIndex = nextIndex + 1;
       end
    else   % check for finished condition in case of fEvery false
      if (t >= tLast)
        % interpolate and exit the while
        d1 = (yNodes(:, 2) - yNodes(:, 1))/(tNodes(2) - tNodes(1));
        d2 = (y - yNodes(:, 2))/(t - tNodes(2));
        d3 = (d2 - d1)/(t - tNodes(1));
        y = yNodes(:, 1) + d1*(tLast - tNodes(1) )+ d3*(tLast - tNodes(1))*(tLast - tNodes(2));
        break;    % we're done, so break out of the main while loop
      end;
    end % if fEvery
    % if we're here, we're not done, so do the updates
    yNodes(:, 1) = yNodes(:, 2);
    yNodes(:, 2) = y;
    tNodes(:, 1) = tNodes(:, 2);
    tNodes(:, 2) = t;
    if (kappa < h*eMin)
      % error is too small so size up h for next round
      h = S*h;
    end;
    else  % error unacceptable, so reject new values
      h = s*h;
      if (h < hMin)
        % stepsize too small, so break out of the while
        gError = true;
        break;
      end;
  end; % end acceptable error
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
