Download<< Back
function [t,x] = ms2(field, t0, T, x0, h)
%AM2 Implements the Milne-Simpson two-step ODE solver
%
% Parameters:
% field -- Right hand side function of ODE system: x'=f(t,x)
% t0 -- Initial time
% T -- End time (T > t0)
% x0 -- Initial value
% h -- Size of time step (h <= T-t0)
%
% Outputs:
% t -- [t0; t-0+h, t0+2*h; ...; t0+i*h; ...]
% x -- Column vector containing numerical solution at each time step
tol = 0.05;
q = 1;
n = ceil((T-t0)/h);
t = t0+h*(0:n).';
x = ones(n+1,length(x0));
[t(1:(q+1)), x(1:(q+1),:)] = rk_classical(field, t0, t0+q*h, x0, h);
f0 = feval(field, t(q+1), x(q+1,:)');
f1 = feval(field, t(q), x(q,:).');
for i=(q+1):n
f_old = feval(field, t(i)+h, x(i,:).');
x(i+1,:) = x(i-1,:).' + h*(f_old+4*f0+f1)/3;
f_new = feval(field, t(i)+h, x(i+1,:).');
while norm(f_old-f_new, 'inf') >= tol
f_old = f_new;
x(i+1,:) = x(i-1,:).' + h*(f_old+4*f0+f1)/3;
f_new = feval(field, t(i)+h, x(i+1,:).');
end
f1 = f0;
f0 = f_new;
end
