function [x, delta, its] = broyden (F, x, tol, max_its)
%
% [x, delta, its] = broyden (F, x0, tol, max_its)
%
% Apply Broyden's method to find the zeros of the vectorial function F.
% x0 is 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. If
% tol is not given, tol is set to [1e-6, 1e-6]. If max_its is not
% given, it is set to 20. delta is an (over)-estimate of the absolute
% error and its the total number of iterations.
if ((nargin <= 3) || isempty (tol))
   tol = [1e-6, 1e-6];
elseif (isscalar (tol))
  tol = [tol, tol];
end
if (nargin <= 4)
  max_its = 20;
end
B = eye (length (x));
Fx = F(x);
delta = -B \ Fx;
its = 1;
while ((norm (delta, inf) > tol(1) * norm (x, inf) + tol(2)) && ...
       (its < max_its))
  x = x + delta;
  y = F(x) - Fx;
  B = B + (y - B * delta) * delta' / (delta' * delta);
  Fx = F(x);
  delta = -B \ Fx;
  its = its + 1;
end
x = x + delta; % extra iteration
delta = norm (delta, inf);
if (delta> tol(1) * norm (x, inf) + tol(2))
  warning ('broyden: tolerance not met')
end
%!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, delta, its] = broyden (F, x0, tol, 8);
%! assert (x, [4.76095822537554e-01; -8.79393408982743e-01], 1e-11)
