% program developed by M. R. Grasselli & D. E. Pelinovsky
% as part of the courseware "Numerical Interpolation &
% Approximation Theory" (2006)

% script_multi_Newton
tol=0.0000001; total = 100; 
% The 'solve' function required the symbolic toolbox which is not
% available on some computer, hence we hardwire the result;
%[x0,y0]=solve('y^2-x=2','x*y=1')
%[sol]=eval([x0,y0])

x0 = [-1;  -1/2+1/2*5^(1/2); -1/2-1/2*5^(1/2)];
y0 = [-1;   1/2*5^(1/2)+1/2;  1/2-1/2*5^(1/2)]; 
sol=[x0,y0];

% iterations of the Newton method, first root
x=[-1.5;-1.5];k = 0;epsilon=1;
while ((epsilon > tol) & (k < total))
    F = [x(1)*x(2)-1;x(2)^2-x(1)-2];
    DF=[x(2) x(1);-1 2*x(2)];
    xnew = x-DF\F;
    epsilon = norm(xnew-x);
    Error = norm(xnew-sol(1,:)');
    x=xnew;
    k = k+1;
fprintf('x = %12.10f, y=%12.10f, k= %2.0f, Error= %12.10f \n', ...
    x(1),x(2),k,Error);
end
fprintf('\n');
% iterations of the Newton method, second root
x=[1;2]; k = 0;epsilon=1;
while ((epsilon > tol) & (k < total))
    F = [x(1)*x(2)-1;x(2)^2-x(1)-2];
    DF=[x(2) x(1);-1 2*x(2)];
    xnew = x-DF\F;
    epsilon = norm(xnew-x);
    Error = norm(xnew-sol(2,:)');
    x=xnew;
    k = k+1;
fprintf('x = %12.10f, y=%12.10f, k= %2.0f, Error= %12.10f \n', ...
    x(1),x(2),k,Error);
end
fprintf('\n');
% iterations of the Newton method, third root
x=[-1.5;-0.5]; k = 0;epsilon=1;
while ((epsilon > tol) & (k < total))
    F = [x(1)*x(2)-1;x(2)^2-x(1)-2];
    DF=[x(2) x(1);-1 2*x(2)];
    xnew = x-DF\F;
    epsilon = norm(xnew-x);
    Error = norm(xnew-sol(3,:)');
    x=xnew;
    k = k+1;
fprintf('x = %12.10f, y=%12.10f, k= %2.0f, Error= %12.10f \n', ...
    x(1),x(2),k,Error);
end
