function x = newtonsys(f,x0)
% NEWTONSYS   Newton's method for a system of equations.
% Input:
%   f        function that outputs values of f(x) and the jacobian
%   x0       initial root approximation
% Output       
%   x        array of approximations (one per column, last is best)

% Operating parameters.
funtol = 100*eps;  xtol = 100*eps;  maxiter = 40;

x = x0(:);  [fn,Jn] = f(x0);
dx = Inf;
n = 1;

while (norm(dx) > xtol) & (norm(fn) > funtol)
  dx = -Jn \ fn;   % Newton step
  x(:,n+1) = x(:,n) + dx;

  n = n+1;
  if n==maxiter
    warning('Maximum number of iterations reached.')
    break
  end

  [fn,Jn] = f(x(:,n));
end