function [x,F,J] = newbat(func,x0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prototype: [x,F,J] = newbat(func,x0)
%
% Purpose:   Just applies Newton's method to the problem F(x)=0,
%            using sufficient descent backtracking linesearch for robustness.
%
% Input:     x0   = initial guess at the minimizer.
%            func = string containing the name of a m-file that
%                   defines the function F(x) and its jacobian J(x),
%                   with prototype: [F,J] = func(x)
%
% Output:    x     = solution of F(x)=0
%            F,J   = F and J at the solution x
%
% Author:  M. Holst (implementation follows P. Gill's notes closely)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tol   = sqrt(eps);    % tolerance for ||F||_2
itmax = 50;           % limit on function evaluations
dpmax = 20;           % limit on damping steps
mu    = 1/4;          % sufficient decrease paramter
delta = 1.0e-10;      % smallest eigenvalue cut-off
iter  = 0;              
alpha = 0;
dmin  = 0;
x     = x0;   
[F,J] = feval(func,x);
normF = norm(F);
fprintf(' iter alpha         ||F||_2\n' )
fprintf(' %3g %9.2e %14.7e\n', iter, alpha, normF);
while ( (normF >= tol) & (iter < itmax) )
   p           = - J\F;
   alpha       = 1;
   xnew        = x + alpha*p;
   [Fnew,Jnew] = feval( func,xnew );
   normFnew    = norm( Fnew );
   j = 1;
   while ( (normFnew > (1-alpha*mu)*normF) & (j <= dpmax) )
      alpha = alpha/2;
      xnew  = x + alpha*p;
      [Fnew,Jnew] = feval(func,xnew);
      normFnew = norm(Fnew);
      j = j + 1;
   end
   x     = xnew;  
   F     = Fnew;        
   J     = Jnew;
   normF = normFnew;
   iter  = iter + 1;
   fprintf(' %3g %9.2e %14.7e\n', iter, alpha, normF);
end
