function  [x,f,k] = Newton_amortiguado(fun,x)
% Método de Newton amortiguado con paso completo
%
 global h
 k=0;  kmax=500;  eps1=1e-9; eps2=1e-12; n=length(x); h=eps1;
 [f,g,H] = fun(x); ng = norm(g,inf); mu = 1; found = ng<=eps1; 
 while ~found
    posdef = 0;
    while  ~posdef
       [R,p] = chol(H+mu*eye(n));   %  Cholesky comprueba si H es definida positiva
       if  p, mu=2*mu; else  posdef = 1; end
    end
    p = -R\(R'\g);  nh = norm(p);    %  Cálculo de paso
    if norm(g,inf)<eps1 || nh <= eps2*(eps2 + norm(x)), found=2;    % Óptimo
    else                             %  Modificar región de confianza
       xnew = x + p;  dL = -p'*g-p'*(H*p)/2;
       [fn,gn,Hn] = fun(xnew);  df = f-fn;  r=df/dL;
       if  r>0
          x = xnew;  f = fn;  g = gn;  H = Hn;
          ng = norm(g,inf);
          mu = mu*max(1/3,1-(2*r-1)^3);
       else
          mu = mu*2;
       end
       k = k + 1;
       if ng <= eps1, found = 1; elseif  k > kmax, found = 3; end
       fprintf('%4.0f %13.8e %13.8e %13.8e %13.8e\n',k,x,f,mu);
    end
 end
end

function [f,g,H]=fun5_2(x)
 f=0.5*x(1)^2*(x(1)^2/6+1)+x(2)*atan(x(2))-0.5*log(x(2)^2+1);
 g=[x(1)^3/3+x(1); atan(x(2))];
 H(1,1)=x(1)^2+1;
 H(1,2)=0;
 H(2,1)=0;
 H(2,2)=1/(1+x(2)^2);
end