
function [xs] = newtonds(x0,f,Tol,MaxIt);
%-----------------------------------------------------------------------
%             -------------                             -------
%             function file                             newtonds.m
%             -------------                             -------
%-----------------------------------------------------------------------
%                        ----------------------------
%                        Zeros Reais de Funcoes Reais
%                        ----------------------------
%                         -------------------------------
%                         Metodo de Newton Para Sistemas
%                          -----------------------------
%
%            Chamada do Procedimento  newtonds
%-----------------------------------------------------
%       xs = newtonsd(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 with damping for nonlinear systems

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

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

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

  for k=1:MaxIt
    if (Defect<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;
      Error = 'Jacobian is singular';
      return;
     else
      dxk = DF\Fk;
    end;
                              
    mt = 1;
    while (mt==1)
      % Compute trial iterate
      xkp1 = xk-theta*dxk;

      %Evaluate function
      Fkp1 = feval(f,xkp1);
      fCalls = fCalls+1;
      mt = (norm(Fkp1,2)>=Defect);

      % Compute a posteriori damping factor
      if (mt==1) 
        theta = 0.5*theta;
       else 
        theta = min(1.0,2.0*theta);
      end;

      if (norm(dxk,2)<1e-12) 
        Error = 'No further change of solution';
        xs = xk;
        return;
      end;
      if (theta<1e-16) 
        Error = 'No further change of solution';
        xs = xk;
	return;
      end;
    end;

    xk = xkp1;
    Fk = Fkp1;
    Defect = norm(Fk,2);
    ItNr = ItNr+1;
  end;

  Error = 'Too many iterations needed';
  xs = xk;
% end of function newtonsd
