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