% Linear multistep method (only k=2)
% Experimental code 

%define the problem
f = @(t,y)  -2*y*t;
ye = @(t) exp(-t.^2)    % Exact solution 
t0 = 0;
tend = 1;
y0 = 1;

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

% Starting values
t1 = t0+h;
y1 = ye(t1);

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

%a = [-0,-1,1];
%b = [-0.5,1.5,0];

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

subplot(2,1,1)
plot(tn,yn,'.-',tn,ye(tn),'--');  
legend('Numerical','Exact')
ylabel('y')
subplot(2,1,2)
plot(tn,ye(tn)-yn);
xlabel('t')
ylabel('Error')

