% Broyden05Test.m
% 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, 04/14/05

clear all; format long; format compact;
maxIt = 200;                % max. number of iterations allowed
tol2 = eps;                 % tolerance in the step size
it = 0;                     % iteration counter                   
Xvec = [0.1; 0.1; -0.1];    % initial guess of root

% Function whose root we are computing
fcn = 'BFp611'; 
% Initial approximation for the inverse of the Jacobian matrix
% Here we use the exact value for faster convergence 
JINV = inv(feval('BFp611Jacobian',Xvec)); 
% In the absence of info on the inverse of the Jacobian matrix,
% one can use the identity matrix. 
JINV = 0.1*eye(length(Xvec));   % better to put in the 0.1 factor since the initial JINV is quite inaccurate for this example!
Vvec = feval(fcn,Xvec); 
Step = -JINV*Vvec;          % step
StepTranspose = transpose(Step);
n2 = StepTranspose*Step;    % Norm squared
Xvec = Xvec+Step;           % Next approximation to the root
it = it+1;
disp(['     loop                               transpose(Xvec)                         |step|']);
while it < maxIt & n2 > tol2
    Wvec = Vvec;
    Vvec  = feval(fcn,Xvec);
    Yvec = Vvec-Wvec;
    Zvec = -JINV*Yvec;
    p = -StepTranspose*Zvec;
    UvecTranspose = StepTranspose*JINV;
    JINV = JINV+(Step+Zvec)*UvecTranspose/p;
    Step = -JINV*Vvec;
    StepTranspose = transpose(Step);
    n2 = StepTranspose*Step;
    Xvec = Xvec+Step;
    disp([it transpose(Xvec) sqrt(n2)]);
    it = it+1;
end


