% newtons1.m
function [xs,Defectf,Error1,fCalls] = newtons1(x0,f,Tol,MaxIt);
%-----------------------------------------------------------------------
%             -------------                             -------
%             function file                             newtons1.m
%             -------------                             -------
%-----------------------------------------------------------------------
%                        ----------------------------
%                        Zeros Reais de Funcoes Reais
%                        ----------------------------
%                         -------------------------------
%                         Metodo de Newton Para Sistemas
%                          -----------------------------
%
%            Chamada do Procedimento  newtons
%-----------------------------------------------------
%       xs = newtons1(x0,f,Tol,MaxIt);
%       Variaveis de Entrada
%---------------------------
%
%       [x0]  :         Vetor com o chute Inicial
%       f :             Funçao
%       tol :           Precisao da Aproximacao do zero da funcao(default=1e-6)
%       MaxiT:          Numero Maximo de Iteraçoes(default=250)
%       Variaveis de Saida
%-------------------------
%
%       xs       :       Valor do zero da funcao
%-------------------------------------------------------------------

if nargin<3
    Tol=1e-6;
end
if nargin<4
    MaxIt=250;
end


% Modified Newton method for nonlinear systems

  h      = 1e-5;       % for approximation of derivative
  n      = max(size(x0));
  Error1 = 'No errors';
  xk     = x0;
  fCalls = 0;

  DF = zeros(n,n);
  Fk = feval(f,xk);    % evaluation of function f like parameter
  fCalls = fCalls+1;
  Defectf = norm(Fk,2);
  ndxk = 1e12;

  for k=1:MaxIt
    if (ndxk<Tol) 
      Error1 = 'No further change of solution';
      xs = xk;
      return;
    end;
    if (Defectf<Tol) 
      xs = xk;
      return;
    end;

    % Evaluate approximate Jacobian
    for j=1:n
      e = zeros(n,1);
      e(j) = 1.0;
      DF(1:n,j) = (1/h)*(feval(f,xk+h*e)-Fk);
    end;
    fCalls = fCalls+n;
       
    % Compute ordinary Newton correction
    if (rcond(DF)>1e12)
      xs = xk;
      Error1 = 'Jacobian is (near) singular';
      return;
     else
      dxk = DF\Fk;
    end;
                              
    xkp1 = xk-dxk;

    %Evaluate function
    Fkp1 = feval(f,xkp1);
    fCalls = fCalls+1;

    xk = xkp1;
    Fk = Fkp1;
    ndxk = norm(dxk,2);
    Defect = norm(Fk,2);
  end;

  Error1 = 'Too many iterations needed';
  xs = xk;
% end of function newtons1
