% Broyden method for the roots of 
% a set of nonlinear equations
% Requires no explicit computation of Jacobian
% See example 5.16 on p.242 of Heath
% This version requires only O(n^2) arithmetic operations/loop.
% K. Ming Leung, 02/21/03

clear all; format long;
profile on;
maxIt = 200; tol2 = eps;
x=[1; 2];
it = 0; err = 1;
hundredK=100000;

tic;
for i = 1:hundredK

% B = feval('fcnJEg5p15',x); A=inv(B);;
% f = feval('fcnEg5p15',x);
% s = -A*f;     % Change in the root

% Much faster (a factor of 2) overall by replacing the 
% above 3 lines with the following line.
% Requires no Jacobian (n^2 function evaluations)
% nor a matrix inverse (O(n^3) arith. ops.)
% Although the number of iterations has to be doubled.
A=eye(length(x)); s=-feval('fcnEg5p15',x); 

n2 = s.'*s;   % Norm squared
x = x+s;      % Next approximation to the root
while it < maxIt & n2 > tol2
    f = feval('fcnEg5p15',x);;
    u = A*f;
    v = s.'*A;
    A = A-u*v/(n2+s.'*u);
    s = -A*f;
    n2 = s.'*s;
    x = x+s;
    it = it+1;
end
end
toc;
profile viewer;

