function [t, Y] = bdf2(N)
a = 0;
b = 1;
h = (b-a) / N;
t = linspace(a, b, N+1);
Y = zeros(1, N+1);
f = @(t, y) -2*t*y^2;
df = @(t, y) -4*t*y;

Y(1) = 1;
% primo passo con Eulero
Y(2) = Y(1) + h*f(t(1), Y(1));

for n = 1:N-1
   %calcolo Y(3),Y(4),...Y(N+1)
   %con BDF
   g=@(y) y - 4/3*Y(n+1) + 1/3*Y(n) ...
        - 2/3*h*f(t(n+2), y);
   dg = @(y) 1 - 2/3*h*df(t(n+2), y);
   % Y(n+2) = newton(g, dg, 5, 1);
   K = 1; % o meglio Y(n+1)
   for k = 1:5
       K = K - g(K) / dg(K);
       g(K)
   end
   Y(n+2) = K;
end



...