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

%*****************************************************************************80
%
%% adapt3 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,1) = t0;
  y(i,:) = y0(:);

  while ( t(i,1) < tstop )
  
    y1 = y(i,:) + dt * ( dydt ( t(i,1), y(i,:) ) )';
 
    th = t(i,1)   + dt / 2.0;
    yh = y(i,:)   + dt / 2.0 * ( dydt ( t(i,1), y(i,:) ) )';
    y2 = yh       + dt / 2.0 * ( dydt ( th, yh ) )';

    e = y2 - y1;

    t(i+1,1) = t(i,1) + dt;
    y(i+1,:) = 2.0 * y2 - y1;
    i = i + 1;

    alpha = sqrt ( tol / norm ( e ) );
    alpha = max ( alpha, 0.3 );
    alpha = min ( alpha, 2.0 );
    dt = dt * 0.9 * alpha;
%
%  If we are going past the endpoint, adjust dt.
%
    dt = min ( dt, tstop + eps - t(i,1) );
%
%  Don't let dt get too small!
%
    dt = max ( dt, dt_min );

  end

  return
end

