function [tvec,yvec] = BDF2System(t0,y0,y1,f,h,N,tol)
% [tvec,yvec] = BDF2System(t0,y0,y1,f,h,N,tol)
% Second-order BDF method
% Inputs
% t0,y0,y1: initial condition
% f:     names of the right-hand side function f(t,y)
% h:     stepsize
% N:     number of steps
% tol:   tolerance for fixed point iteration
% Outputs
% tvec: vector of t values
% yvec: vector of corresponding y values

d = length(y0);
yvec = zeros(N+1,d);
tvec = linspace(t0,t0+N*h,N+1)';
yvec(1,:) = y0';
yvec(2,:) = y1';
w = yvec(2,:);
for n = 1:N-1
   converged = 1;   k = 1;
   while  (converged > tol) & (k < 10)     % fixed point iteration to solve implicit problem
       fvalue = f(tvec(n+2),w);
       w = (4*yvec(n+1,:)-yvec(n,:)+2*h*fvalue')/3;
       if k > 1
           converged = abs(wold-w);
       end
       wold = w;
       k = k+1;
   end    
   yvec(n+2,:) = w;
end