function [xx, it] = broyden(x, f, n, tol)
fx = zeros(n,1);
it = 0;
xx = x;
B = eye(n);
fx = feval(f, xx);
while norm(fx) > tol
   it = it + 1;
   h = - B*fx;
   lambda = 1;
   xx = xx + lambda*h;
   fprintf('Iteration %d: approx root = [%f,%f]\n', it, xx(1), xx(2));
	oldfx = fx;
   fx = feval(f, xx);
   z = fx - oldfx;
   oldB = B;
   ozh = oldB*z - h;
   hB = h'*oldB;
   for i = 1:n
      for j = 1:n
         M(i,j) = ozh(i)*hB(j);
      end;
   end;
   B = oldB - M./(h'*oldB*z)
end;
