function Levmar_99
% Mínimos cuadrados con Levenberg-Marquardt
%
 x=[10;1;1]; mu=0.01; dnor=1; n=3;
 tol=sqrt(eps)*norm(x,inf); update=1; k=0; i=0;
 while dnor>tol
    if update==1
       f=fx(x);
       J=derf(x);
       jtj=J'*J;
       res=norm(f)^2;
    end
    a=jtj+mu*eye(n);
    s=a\(J'*f);
    b=x-s;
    f1=fx(b);
    res1=norm(f1)^2;
    if res1<res
       x=b;
       f=f1;
       dnor=norm(s,inf)/norm(x,inf); k=k+1;
       fprintf('%3.0f %15.10e %15.10e %15.10e %15.10e %15.10e %15.10e\n',...
          k,x,res1,mu,dnor);
       mu=mu/10;
       update=1;
    else
       mu=mu*10;
       update=0; i=1+1;
    end
 end
end

function f=fx(x)
 f     = zeros(12,1);
 f(1)  = x(1)/(1+x(2)*exp(x(3)))-5.308;
 f(2)  = x(1)/(1+x(2)*exp(2*x(3)))-7.24;
 f(3)  = x(1)/(1+x(2)*exp(3*x(3)))-9.638;
 f(4)  = x(1)/(1+x(2)*exp(4*x(3)))-12.866;
 f(5)  = x(1)/(1+x(2)*exp(5*x(3)))-17.069;
 f(6)  = x(1)/(1+x(2)*exp(6*x(3)))-23.192;
 f(7)  = x(1)/(1+x(2)*exp(7*x(3)))-31.443;
 f(8)  = x(1)/(1+x(2)*exp(8*x(3)))-38.558;
 f(9)  = x(1)/(1+x(2)*exp(9*x(3)))-50.156;
 f(10) = x(1)/(1+x(2)*exp(10*x(3)))-62.948;
 f(11) = x(1)/(1+x(2)*exp(11*x(3)))-75.995;
 f(12) = x(1)/(1+x(2)*exp(12*x(3)))-91.972;
end

function J=derf(x)
J       = zeros(12,3);
J(1,1)  = 1/(1+x(2)*exp(x(3)));                             
J(1,2)  = -x(1)*exp(x(3))/(1+x(2)*exp(x(3)))^2;
J(1,3)  = -x(1)*x(2)*exp(x(3))/(1+x(2)*exp(x(3)))^2;           
J(2,1)  = 1/(1+x(2)*exp(2*x(3)));          
J(2,2)  = -x(1)*exp(2*x(3))/(1+x(2)*exp(2*x(3)))^2;            
J(2,3)  = -x(1)*x(2)*exp(2*x(3))*2/(1+x(2)*exp(2*x(3)))^2;     
J(3,1)  = 1/(1+x(2)*exp(3*x(3)));    
J(3,2)  = -x(1)*exp(3*x(3))/(1+x(2)*exp(3*x(3)))^2;            
J(3,3)  = -x(1)*x(2)*exp(3*x(3))*3/(1+x(2)*exp(3*x(3)))^2;     
J(4,1)  = 1/(1+x(2)*exp(4*x(3)));    
J(4,2)  = -x(1)*exp(4*x(3))/(1+x(2)*exp(4*x(3)))^2;            
J(4,3)  = -x(1)*x(2)*exp(4*x(3))*4/(1+x(2)*exp(4*x(3)))^2;     
J(5,1)  = 1/(1+x(2)*exp(5*x(3)));    
J(5,2)  = -x(1)*exp(5*x(3))/(1+x(2)*exp(5*x(3)))^2;            
J(5,3)  = -x(1)*x(2)*exp(5*x(3))*5/(1+x(2)*exp(5*x(3)))^2;     
J(6,1)  = 1/(1+x(2)*exp(6*x(3)));    
J(6,2)  = -x(1)*exp(6*x(3))/(1+x(2)*exp(6*x(3)))^2;            
J(6,3)  = -x(1)*x(2)*exp(6*x(3))*6/(1+x(2)*exp(6*x(3)))^2;     
J(7,1)  = 1/(1+x(2)*exp(7*x(3)));    
J(7,2)  = -x(1)*exp(7*x(3))/(1+x(2)*exp(7*x(3)))^2;            
J(7,3)  = -x(1)*x(2)*exp(7*x(3))*7/(1+x(2)*exp(7*x(3)))^2;     
J(8,1)  = 1/(1+x(2)*exp(8*x(3)));    
J(8,2)  = -x(1)*exp(8*x(3))/(1+x(2)*exp(8*x(3)))^2;            
J(8,3)  = -x(1)*x(2)*exp(8*x(3))*8/(1+x(2)*exp(8*x(3)))^2;     
J(9,1)  = 1/(1+x(2)*exp(9*x(3)));                                
J(9,2)  = -x(1)*exp(9*x(3))/(1+x(2)*exp(9*x(3)))^2;            
J(9,3)  = -x(1)*x(2)*exp(9*x(3))*9/(1+x(2)*exp(9*x(3)))^2;     
J(10,1) = 1/(1+x(2)*exp(10*x(3)));                               
J(10,2) = -x(1)*exp(10*x(3))/(1+x(2)*exp(10*x(3)))^2;          
J(10,3) = -x(1)*x(2)*exp(10*x(3))*10/(1+x(2)*exp(10*x(3)))^2;  
J(11,1) = 1/(1+x(2)*exp(11*x(3))); 
J(11,2) = -x(1)*exp(11*x(3))/(1+x(2)*exp(11*x(3)))^2;          
J(11,3) = -x(1)*x(2)*exp(11*x(3))*11/(1+x(2)*exp(11*x(3)))^2;  
J(12,1) = 1/(1+x(2)*exp(12*x(3))); 
J(12,2) = -x(1)*exp(12*x(3))/(1+x(2)*exp(12*x(3)))^2;          
J(12,3) = -x(1)*x(2)*exp(12*x(3))*12/(1+x(2)*exp(12*x(3)))^2;
end
