
function [xs,ItNr] = newtons(x0,f,Tolx,Tolf,MaxIt);
%-----------------------------------------------------------------------
%             -------------                             -------
%             function file                             newtons.m
%             -------------                             -------
%-----------------------------------------------------------------------
%                        ----------------------------
%                        Zeros Reais de Funcoes Reais
%                        ----------------------------
%                         -------------------------------
%                         Metodo de Newton Para Sistemas
%                          -----------------------------
%
%            Chamada do Procedimento  newtons
%-----------------------------------------------------
%       xs = newtons(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
%-------------------------------------------------------------------

% Modified Newton method for nonlinear systems
if nargin<3
    Tolx=1e-6;
end
if nargin<4
    Tolf=1e-6;
end
if nargin<5
    MaxIt=250;
end

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

  DF = zeros(n,n);
  Fk = feval(f,xk);    % evaluation of function f like parameter
  Defectf = norm(Fk,2);
  Defectx = 0;

  for k=1:MaxIt
    if (Defectf<Tolf) 
      Error1 = 'Small norm value of function';
      xs = xk;
      return;
    end;
    % Evaluate Jacobian without its control of regularity
    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;
    % Compute ordinary Newton correction
    dxk = DF\Fk;
    % Compute trial iterate.
    xkp1 = xk-dxk;
    ItNr = ItNr+1;
    Defectx = abs(dxk); 
    % Evaluate function
    Fkp1 = feval(f,xkp1);
    if (norm(dxk,2)<Tolx) 
      Error1 = 'No further change of solution';
      xs = xk;
      return;
    end;

    xk = xkp1;
    Fk = Fkp1;
    Defectf = norm(Fk,2);
  end;

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