% Linear multistep method

%define the problem
f = inline('-2*t*y','t','y');
t0 = 0;
tend = 1;
y0 = 1;

Nstep = 10;
h = (tend-t0)/Nstep;

y1 = exp(-h^2);
t1 = t0+h;

% Method (explicit, k=2);
a = [-1,0,1];      % [alpha0,alpah1,alpha2]
b = [0,2,0];       % [beta0,beta1,beta2]

yn = [y0,y1];
tn = [t0,t1];
fn = [f(t0,y0),f(t1,y1)];
for n=2:Nstep 
    y = -a(2)*yn(:,n)-a(1)*yn(:,n-1) ...
      + h*(b(2)*fn(:,n)+b(1)*fn(:,n-1));
    t = tn(n)+h;  
    yn = [yn,y];
    tn = [tn,t];
    fn = [fn,f(t,y)];
end

ye = exp(-tn.^2);
plot(tn,yn,'.-',tn,ye,'--');  

