function [t,y] = fadapt(FunFcn,tspan,y0,abstol,reltol,hmax,hmin)
% This function uses Heun's method and Euler's method to adaptively produce
% an approximation to the solution of the system of differential equations
% y' = f(t,y),  y(t0) = y0
% In the calling sequence
% [t,y] = adapt('f',tspan,y0,abstol,reltol,hmax,hmin)
% 'f' is a string containing the name of matlab m file which contains the
% function f(t,y).  This function takes as inputs a time t and a column vector
% y and returns the column vector evaluating f(t,y).
% tspan = [t0, tfinal] is a vector of the initial and final times
% The initial condition y0 must be a column vector
% Each row in the solution array y corresponds to a time returned in column
% vector t, i.e., y(i,j) gives an approximation to y_j(t(i)).
% abstol and reltol are absolute and relative error tolerances.
% The goal is to produce an approximation y such that the approximation
% est to the local error per unit step satisfies 
% est_j <= max(rtol*abs(y_j),atol) for each component j of the solution.
%
% Initialize: Check to see if abstol, reltol, hmax, and hmin are provided
% by the user.  If not, use default values.  Note that
% isempty(x) returns 1 if x is an empty matrix and 0 otherwise
% In the case a non-zero element is returned, the following if statement
% is executed.
%
if isempty(abstol) atol= 1e-4
else atol=abstol;
end
if isempty(reltol) rtol = 1e-3
else rtol=reltol;
end
t0=tspan(1);
tfinal=tspan(2);
if isempty(hmax) hmax = (tfinal-t0)/10;
end
if isempty(hmin) hmin  = 16*eps*abs(tfinal-t0);
% eps is a permanent constant ~ 2.22e-16
end
h= sqrt(hmax*hmin);   % This is just an arbitrary way of getting a initial h
i=0;
yold = y0;
told = t0;
tnew = told + h;
% General Step
while tnew <= tfinal
% Get approximations by Heun and Euler
k1 = feval(FunFcn, told, yold);    % Note k1 is a column vector
yeuler = yold + h*k1;              % yeuler is a column vector
yheun = yold + .5*h*k1 + .5*h*feval(FunFcn, tnew, yeuler);  % column vector
reject = 0;
% Decide whether to accept or reject this step
% length(v) produces the number of components in the vector v
for j=1:length(yheun)
% compute est(j), the jth component of an estimate of the local error/unit step
     est(j) = abs(yheun(j)-yeuler(j))/h;
     tol(j) = max(rtol*abs(yheun(j)),atol);
     gam(j) = tol(j)/est(j);
end
% Note that the above loop could be replaced by vector equations
% Check whether all components of est are <= tol; if not set reject =1
if max(est-tol) > 0 reject =1;   % max(x) = maximum component of vector x
end
% Compute gamma = min_j[tol(j)/est(j)];
gamma = min(gam);      % min(x) = minimum component of vector x
if reject == 0 
   i = i+1;
   tr(i) = tnew;
   y(i,:) = yeuler';   % the symbol ' is for transpose.  Here it
%                       converts the column vector yeuler to a row vector.
%                       The notation y(i,:) says to place the row vector
%                       yeuler' in the ith row of y.
   told = tnew;
   yold = yeuler;
   ht = min(h*0.8*gamma,hmax);
   h =  min(ht,5*h);  % 0.8 is an arbitrary fudge factor used
% to help avoid rejection of the next step. Also do not allow step size
% greater than hmax and do not increase step size by more than a factor of 5
else
   ht = max(h*0.8*gamma,hmin);
   h =  max(ht,h/10);  % Do not allow step size smaller than hmin
%  and do not decrease step size by more than a factor of 10
end
   tnew = told + h;
end
t = tr';  % The symbol ' converts the row vector tr to a column vector t

%% Some items to add
%% 1. Once t > tfinal, recompute last value of y so that it corresponds
%% to tfinal
%% 2. Plot graph of solution as it is being computed to see adaptivity
%% 3. Use a better method for initializing h (see other matlab codes)