Download<< Back
function [t,x] = ode12(field, t0, T, x0)
%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
%
% 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
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
tol = 0.001;
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
