Personal Website:
Scott Congreve

Download<< Back

function [t,x] = bdf2(field, t0, T, x0, h)
%BDF2 Implements the BDF 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.001;
m = 2;
n = ceil((T-t0)/h);

t = t0+h*(0:n).';
x = ones(n+1,length(x0));
[to,xo] = rk_classical(field, t0, t0+(m-1)*h, x0, h);
t(1:m) = to(1:m);
x(1:m,:) = xo(1:m,:);

for i=m:n
    f_old = feval(field, t(i)+h, x(i,:).');
    x(i+1,:) = (4*x(i,:).' - x(i-1,:).' + 2*h*f_old)/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,:) = (4*x(i,:).' - x(i-1,:).' + 2*h*f_old)/3;
        f_new = feval(field, t(i)+h, x(i+1,:).');
    end
end