function  [x,i,f,nfeval] = quasi_newton_1(fun,x,metodo)   
% Método cuasi Newton para minimizar f(x).
%
 rho = 0.01; beta = 0.1;                % Parámetros de line search
 [f,g] = fun(x); H=eye(length(x)); eps1=1e-5; eps2=1e-8;
 maxit = 1000; nfeval = 1;
 for i=1:maxit
    xp = x;  gp = g;
    p = -H*g;  nh = norm(p);            % Nueva dirección: cuasi Newton
    if norm(g,inf)<eps1 || nh <= eps2*(eps2 + norm(x)), break, end
    alpha = 1;
    for k=1:10                          % Backtracking amplitud de paso
       xnew = x+alpha*p; 
       fxnew = fun(xnew); nfeval = nfeval+1;
       if fxnew < f+alpha*rho*g'*p, break
       else alpha=alpha*beta; 
       end
    end
    x = x + alpha*p; s = x - xp; ns = norm(s);
    [f,g] = fun(x); 
    y = g - gp;  ys = y'*s;  nfeval = nfeval+1;
    if  ys > sqrt(eps)*ns*norm(y)                       % Recalcular  H
        v = H*y;  yv = y'*v;
        if  metodo == 1                                 % DFP
            H = H - (v/yv)*v' + (s/ys)*s';
        else                                            % BFGS
            H = H + (1+yv/ys)*(s*s')/ys - (s*v'+v*s')/ys;    
        end
    end  
 end
end