function [x, y] = ODE_Euler_implicit(f,a,b,y0)
n = 1000;
N = 1000;
x = linspace(a,b,n);
delta = x(2) - x(1);
y = zeros(1,n);
y(1) = y0;
for i = 2 : n
    % Newton's step
    y(i) = y(i-1);
    for j = 1 : N
        d = (f(x(i),y(i) + 0.001) - f(x(i),y(i))) / 0.001 * delta - 1;
        y(i) = y(i) - (f(x(i),y(i)) * delta - y(i) + y(i-1)) / d; 
    end
end
end