function x = newton (f, Jf, x0, tol)
%
% x = newton(f,Jf,x0,tol)
%
% Apply Newton's method to find the zeros of the vectorial function f.
% Jf is its Jacobian, x0 an initial guess and tol the vector of
% tolerances (tol(1) is the relative tolerance, tol(2) the absolute
% tolerance). If tol is a scalar, both tolerances are set to tol.
% In GNU Octave you can type
%
% octave:1> test newton
%
% and the code written as a comment at the bottom of the function will
% be executed.
%
if (isscalar (tol))
  tol = [tol, tol];
end
x = x0;
delta = -Jf(x) \ f(x);
while (norm (delta, inf) > tol(1) * norm (x, inf) + tol(2))
    x = x + delta;
    delta = -Jf(x) \ f(x);
end
x = x + delta; % extra iteration
%!test
%! f = @(x) [x(1) ^ 2 + x(2) ^ 2 - 1; sin(pi * x(1) / 2) + x(2) ^ 3];
%! Jf = @(x) [2 * x(1), 2 * x(2); pi / 2 * cos(pi * x(1) / 2), 3 * x(2) ^ 2];
%! x0 = [1; 1];
%! tol = [1e-8, 1e-8];
%! x = newton (f, Jf, x0, tol);
%! assert (x, [4.76095822537554e-01; -8.79393408982743e-01], 1e-11)
