function [ t, y ] = adaptive3 ( dydt, tspan, y0, dt, tol )

%*****************************************************************************80
%
%% rk1_adaptive_double uses an adaptive stepsize.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 October 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    function handle dydt: defines the derivative function of the form
%      value = dydt ( t, y )
%
%    real tspan(2): contains the initial and final times.
%
%    real y0(m): contains the initial condition.
%
%    real dt: the size of the first step.
%
%    real tol: an error tolerance.
%
%  Output:
%
%    real t(n+1): the times.
%
%    real y(n+,m): the estimated solution values.
%
  t0 = tspan(1);
  tstop = tspan(2);
  dt_min = 1.0E-6;

  i = 1;
  t(i) = t0;
  y(i) = y0;

  while t(i) < tstop

    while ( true )
  
      y1 = y(i) + dt * dydt ( t(i), y(i) );
 
      th = t(i) + dt / 2;
      yh = y(i) + dt / 2 * dydt ( t(i), y(i) );
      y2 = yh   + dt / 2 * dydt ( th, yh );

      e = y2 - y1;
 
      if ( abs ( e ) / dt < tol / ( tstop - t0 ) )
        break;
      elseif ( dt < dt_min )

        error ( 'Stepsize %d is smaller than %g\n', dt, dt_min );
      else
        dt = dt / 2.0;
      end
    
    end
  
    t(i+1) = t(i) + dt;
    y(i+1) = 2.0 * y2 - y1;
    i = i + 1;

    dt = dt * 0.9 * min ( max ( sqrt(0.5*tol/e), 0.3 ), 2.0 );

    if ( 4.0 * abs ( e ) / ( 2.0 * dt ) < tol / ( tstop - t0 ) )
%     dt = 2.0 * dt;
    end
    
  end

  return
end

