function [t,y] = euleradapt( f, tspan, y0)
%EULERADAPT Apply adaptive Eulers's method to solve y' = f(y,t), y(a) = y0,
%  on the interval tspan(1) <= t <= tspan(2).
%
%  This uses Euler's method to step, and the modified Euler's method to
%  compute an error bound and adapt the step size.
%
%  Example usage:
%    f = @(t,y) t.*y.^2;
%    [t,y] = euleradapt(f,[0,1],1,20);
%    plot(t,y);
%
    a = tspan(1);
    b = tspan(2);
    
    tolerance = 1.e-3;
    h = (b-a)*tolerance;  % initial step size
    
    j = 1;
    t(1) = a;
    y(1) = y0;
    
    while t(j) < b
        m1 = f(t(j), y(j));  % slope at current point
        m2 = f(t(j)+h, y(j)+h*m1);  % slope at next point (approx)
        
        % Euler's approximation
        w1 = y(j) + h*m1;
        
        % Modified Euler's approximation
        w2 = y(j) + h*(m1+m2)/2;
        
        % Estimate local truncation error
        err = abs( w1 - w2 )/h;
        if err > tolerance
            % Too much error. Repeat with smaller step size
            h = h/2;
        else
            % Error acceptable. Take a step.
            y(j+1) = w1;
            t(j+1) = t(j) + h;
            j=j+1;
            
            if err < tolerance/4
                % Could be taking larger steps
                % For Euler's method, doubling h approximately doubles the
                % truncation error, resulting in err < tolerance/2
                h = 2*h;
            end
        end
    end % while loop
end