function [x, it] = broyden(x0, f)
% BROYDEN Broyden's method for solving a system of
%         n non-linear equations in n variables.
% Input:
%   f        function that computes residual (n-vector)
%   x0       initial root approximation (n-vector)
%
% Output       
%   x        array of approximations (one per column, last is best)
%   it       number of iterations
%
    it = 1; 
    x = x0(:);  
    y = f(x0);
    dx = Inf;
    
    % Operating parameters.
    funtol = 1000*eps;  
    xtol = 1000*eps;  
    maxiter = 60;

    % Set initial Br
    n = length(x0);
    Br = eye(n);
    
    while (norm(dx) > xtol) && (norm(y) > funtol) && (it < maxiter)
      dx = -Br*y;
      x(:, it+1) = x(:, it) + dx;
      oldy = y; 
      y = f(x(:,it+1));
      
      % Update approximation to inverse Jacobian using Broydens formula
      dy = y - oldy; 
      Br = Br - (Br*dy - dx)*(dx'*Br)/(dx'*Br*dy);
      it = it + 1;
    end
end
