%solve F(x) = 0 with Newton
% where F(x) = [x1^2 + 1/4(x2-a)^2 -1; x1^2+y1^2 - 1];
% for different values of a.
%
%
a   =[0.5;1;1.5;3;4];
tol = 1e-10;
nit = 50;
warning off;
for k = 1:length(a);
    x = [1;1];
%
% no convergence with x=[0;0];
%
    flg = 1;
    for t=1:nit
        Jx = [2*x(1) 1/2*(x(2)-a(k));2*x(1) 2*x(2)];
        Fx = [x(1)^2 + 1/4*(x(2)-a(k))^2-1;x(1)^2+x(2)^2-1];
        if (norm(Fx)<tol)
            disp(['converged, a = ',num2str(a(k)), ' in ', num2str(t), ' iterations'])
            disp(x);
            flg = 0;
            break;
        end
        x  = x - Jx\Fx;
     end
     if (flg == 1)
         disp(['no convergence, a = ',num2str(a(k)), ' in ', num2str(nit), ' iterations'])
     end
end