function [x1, err] = twostep(ODE, x0, t0, h, stepper, order)
%
% Note that err is a vector (same size as phase space)
%

% Estimate x1 by one full step.
x1a = stepper(ODE, x0, t0, h);

% Estimate x1 by two half steps.
% [Note that ODE is called at least once more than is necessary].
thalf = t0 + h/2;
xhalf = stepper(ODE, x0, t0, h/2);
x1b = stepper(ODE, xhalf, thalf, h/2);

% Estimate truncation error, and use it to correct x1.
err = (x1b - x1a)/(2^order - 1); % This is a vector if x is.
x1 = x1b + err; % Order of accuracy is 1+order.
