function [t,x] = ode12(field, t0, T, x0, tol)
%ODE12 Implements adaptive timestepping with p=1 (Euler) and p+1=2 (Heun)
%
% Parameters:
%  field  -- Right hand side function of ODE system: x'=f(t,x)
%  t0     -- Initial time
%  T      -- End time (T > t0)
%  x0     -- Initial value
%  tol    -- Tolerance to use for the adaptive timestepping (default=0.001)
%
% Outputs:
%  t  -- [t0; t-0+h1, t0+2*h2; ...; t0+i*hi; ...]
%  x  -- Matrix containing numerical solution, with each row the value of x
%        at each time step

arguments
    field; t0; T; x0;
    tol = 0.001;
end

t = [t0];
x = zeros(1, length(x0));
x(1,:) = x0;
hmax = (T - t0)/16; % change from ode12_2
h = hmax/8;  % change from ode12_2
power = 1/2; % = 1/(p+1)

while t(end) <= T
    k1 = feval(field, t(end), x(end,:).');
    k2 = feval(field, t(end)+h, x(end,:).'+h*k1);
    % Euler: x(i,:).' + h*k1;
    % Heun:  x(i,:).' + h*(k1+k2)/2;
    % => Heun-Euler = h*(k2-k1)/2;
    delta = norm(h*(k2-k1)/2, 'inf');
    if delta <= tol
        % delta less than tolerance, so next timepoint and solution
        % computed by Heun       
        t = [t;t(end)+h];
        x = [x; x(end,:)+h*(k1.'+k2.')/2];
    end
    if delta ~= 0.0
        h = min(hmax, 0.9*h*(tol/delta)^power);
    end
end
