function [x, res] = dampnewton(x,F,DF,abstol,reltol)
% Damped Newton method for non-linear system of equations
% x passes the initial guess, F and DF handles of type @(x) to the function
% F and its Jacobian, abstol and reltol specify absolute and relative
% tolerances for termination. Returns convergence history beside
% approximate solution.
LMIN = 1e-3; % minimal damping factor
% Compute first Newton correction
[L,U] = lu(DF(x)); s = U\(L\F(x)); 
% first update, initially no damping
xn = x-s; lambda = 1; 
% first simplified Newton correction
f = F(xn); st = U\(L\f); stn = norm(st);
res = stn; % matrix for recording convergence history
while ((stn > reltol*norm(xn)) && (stn > abstol))
  % natural monotonicity test
  while (norm(st) > (1-lambda/2)*norm(s)) 
    % if not passed reduce damping factor by a factor of 2
    lambda = lambda/2;
    % If damping becomes too strong, bail out
    if (lambda < LMIN), warning('No convergence'); return; end
    % New tentative iterate and associated simplified Newton correction  
    xn = x-lambda*s; f = F(xn); 
    st =  U\(L\f);
  end
  % Step accepted, update approximate solution
  x = xn; [L,U] = lu(DF(x)); s = U\(L\f);
  % Boldly increase damping factor
  lambda = min(2*lambda,1);
  % Compute next tentative iterate and associated simplified Newton 
  % correction 
  xn = x-lambda*s; f = F(xn); st = U\(L\f);
  stn = norm(st);
  res = [res; stn]; % record norms of simplified Newton corrections
end
x = xn;
