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;