% Newton's method and Broyden's method for nonlinear systems of equations
TOL=1.e-8; N0=20;
n=3;   % for an example of 3x3 system.

% You can choose one of the two methods: Newton or Broyden
if exist('method') == 0
    method = 'newton';
end

x = zeros(n,1);

switch lower(method)
    
    case 'newton'
    k=1;
    while k <= N0 
        pause
        [f,J] = fvec(x,1);    %In fvec with the flag set to 1, the Jacobian
                              %is computed
        y = - J\f;            %correction
        x = x + y;            %updated approximation
        fprintf('approximation after %d iterations\n',k);
        disp(x);
        if norm(y,inf) < TOL
            fprintf('approximation satisfactory!\n');
            return;
        end
        k = k+1;
    end 
    fprintf('The method failed after %d iterations.\n',N0);
    
    case 'broyden'
    [v,J] = fvec(x,1);        %for the first step, the Jacobian is used.
    Ainv = inv(J);
    s = -Ainv*v;              %correction in the first iteration
    x = x + s;                %first updated solution
    k = 2;
    while k <= N0
        pause
        w = v;                %storing the previous f values
        v = fvec(x,0);        %new f values, but no Jacobian
        y = v - w;            %change in f values
        z = -Ainv*y;          %matrix-vector multiplication
        p = -s'*z;            %inner product
        u = Ainv'*s;          %another matrix-vector multiplication, notice the transpose
        Ainv = Ainv + (s+z)*u'/p;  %correction to A inverse involving a rank-1 matrix
        s = -Ainv*v;          %correction in approximation
        x = x + s;            %updated approximation
        fprintf('approximation after %d iterations\n',k);
        disp(x);
        if norm(s,inf) < TOL
            fprintf('approximation satisfactory!\n');
            return;
        end
        k = k+1;
    end               
    otherwise
        disp('method must be either Newton or Broyden.');

end
return