function retval = Trap(fName, t0, y0, hStep, nSteps)
% usage: retval = Trap(fName, vt0, 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 an Euler
% method with Heun (Ralston) method as a starter, followed by
% a set number of fixed point iterations on the trapezoidal method. Default
% for the optional parameter theta is 1/2, ie, trapezoidal method.
% 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: 2/27/16

% Octave and Matlab compatible file
% local variables:
% t: current time value
% jj: index variable
% iterLimit: limit on iterations (must be >= 1 and is highly problem dependent)
% fEvery: flag to return initial and every subsequent profile in a matrix

global gError gxCount; 

% initialization code
gError = false;
gxCount = 0;
% initialization code
fEvery = false;  % set to true for full profile
t = t0;
iterLimit = 2;
newVal= y0;
if (fEvery)
  retval = zeros(length(y0), nSteps + 1);
end
% main loop
for jj = 1:nSteps
  % avoid extra evaluations
  fcn0 = feval(fName, t, y0);
  % do the Euler step
  y1 = y0 + hStep*fcn0;
  % do the trapezoidal step with fixed point iteration
  for ii = 1:iterLimit
    fcn1 = feval(fName, t + hStep, y1);
    y1 = y0 + 0.5*hStep*(fcn0 + fcn1);
  end
  % set up for next iteration
  y0 = y1;
  t = t + hStep;
  if (fEvery)
    retval(:, jj + 1) = y1;
  end  
end;
if (~fEvery)
  retval = y1;
end

end
