function [Y,T] = vareuler(y0,tau,Tend,fun,param,tol);
%
%  [Y,T] = vareuler(y0,tau,Tend,fun,param,tol);
%  
%  Integrates the function 'yprime = fun (y,t,param)'
%  using Forward Euler.
%
%  y0 is the initial condition
%  tau is an initial stepsize guess
%  Tend is the length of the integration interval
%  tol is an absolute tolerance
%
%  The output is stored columnwise in the array Y
%  The array T contains the temporal grid
%

Y = y0;  %% output array
y = y0;  %% solution at current time level
T = 0;   %% output array
t = 0;   %% current time level

facmin = 0.5;
facmax = 2;
fac = 0.8;


while t<Tend,  %% Note:  we do not know in advance how many steps are needed
   remesh = 1; %% Assume remeshing is necessary (1==true, 0==false) 
   
   while remesh
      
      %% Euler step
      F = feval(fun,y,t,param);
      y1 = y + tau * F;
      
      %% Richardson extrapolation
      ya = y + tau/2 * F;
      yb = ya + tau/2 * feval(fun,ya,t+tau/2,param);  
      y2 = 2*yb-y1; 
      
      %% Error estimate:
      err = max(abs(y2-y1))/tol;
   
      if err <= 1
         remesh = 0;  %% remeshing unnecessary, drop out of the loop
         t = t + tau; %% increment time
      end
      
      %% Estimate next timestep
      tau = min(Tend-t,tau * min(facmax,max(facmin,fac/sqrt(err))));
   end
    
   y = y1;
   
   Y(:,end+1) = y;
   T(end+1) = t;
   
end



