function [tt, yy] = ode324(f, tspan, y0, tol)
% ODE324 - ODE IVP sovler by Dr Hale and TW324 class 2017
% INPUTS:  A function handle f = @(t,y) ... such that dy/dt = f(t,y)
%          The times span tspan(1) < t <tspane(end)
%          The initial condition y(tspan(1)) - y0
%          The required accuracy, tol.
% OUTPUTS: The selected time steps tt, and solution values yy.
%          (Note, if no outputs requested, then a plot is shown instead).

% Initial step size
h = sqrt(tol);      
% Initialise solver
t    = tspan(1);
tend = tspan(2);
yn   = y0;
flag = 0;
% Storage
tt = t;
yy = yn;

%% Loop through time
while ( t < tend )
    
    % RK stages
    s1 =                            ;
    s2 =                            ;
    s3 =                            ;
    p  =                            ;
    yn1 = yn +                      ;
    % Error estimate
    err =                           ;
       
    % Step size adjust:
    hnew = 0.8*(tol*min(abs(yn1)./err))^(1/(p+1))*h;
    hnew = min(hnew, 2*h);          % Don't let h grow too fast.
    
    if ( norm(err, inf) > tol  )    % Error is too large:
        if ( ~flag )                % Try the new time step,
            flag = 1;               % Set flag to remember failure.
        else                        % Time step failed. 
            hnew = h/2;             % Take half of old one,
        end
    else                            % Error is small enough!
        % Store solution and move to next time:
        yn = yn1;
        t = t + h;
        yy = [yy, yn1];          
        tt = [tt, t];
        flag = 0;                   % Successful step. Set flag to happy.
    end   
   
    % Choose next step (ensure we finish at the end of the interval!)
    h = min(hnew, tend - t);           

end

if ( nargout == 0 )
    subplot(2,1,1)
    plot(tt, yy(1,:), '.'), shg
    title('First component of solution')
    subplot(2,1,2)
    semilogy(tt(2:end), abs(diff(tt)))
    title(['No. of time steps = ' int2str(length(tt))])
    tt = [];
    yy = [];
end

end
