function [x1,f1,niter] = bfgs(f,x0,h,tolgrad,toldiff,tolline,nder,max_b_iter,max_s_iter)
% f = function to be minimized
% x = starting point as a column vector.
% h = initial search increment in the line search.
% nder = type of derivative: [1] 1st forward [2] 2nd forward [3] central
% tolgrad = error tolerance on the gradient.
% toldiff = error tolerance on the var diff.
% tolline = error tolerance on the minimum of f(x) in the line search.
% max_b_iter = maximum number of iterations in the bracketing phase.
% max_s_iter = maximum number of iterations in the sectioning phase.

if nargin < 9; max_s_iter = 10000; end
if nargin < 8; max_b_iter = 10000; end
if nargin < 7; nder = 3; end
if nargin < 6; tolgrad = 1.0e-8; end
if nargin < 5; toldiff = 1.0e-6; end
if nargin < 4; tolgrad = 1.0e-4; end
if nargin < 3; h = 0.01; end

n = size(x0,1);
I = eye(n);
H0 = I;
[~,g0] = gradx(f,x0,nder);

niter = 0;
s = toldiff;

while norm(s) >= toldiff  
    
        if (g0'*g0) <= tolgrad
            return
        end

        niter = niter + 1;
        p = -H0*g0;
        
        % Make sure the search direction is real.
        % p = real(p);
                    
        % Line search.
        [s,~] = line_search(@(s) f(x0+s*p),0.0,-Inf,Inf,h,tolline,...
                max_b_iter,max_s_iter);
        x1 = x0+s*p;
                
        s = x1 - x0;
        [f1,g1] = gradx(f,x1,nder);
        
        q = g1 - g0;
        rho = 1/(q'*s);
        H1 = (I - rho*s*q')*H0*(I-rho*q*s') + rho*(s*s');
        
        x0 = x1;g0 = g1;H0 = H1; 

end

end
