% Linear multistep method

%define the problem
f = @(t,y) -2*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 = [0,-1,1];           % [alpha0,alpah1,alpha2]
b = [-0.5,1.5,0];       % [beta0,beta1,beta2]

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

plot(tsol,ysol,'o-',tsol,exp(-tsol.^2),'--');  

