function [list_iter,list_steps] = BFGS(f_func, g_func, init_pt, H0, N_iter, epsilon)
list_iter = zeros(length(init_pt),N_iter+1);
list_iter(:,1) = init_pt;
list_steps = zeros(N_iter,1);

x = init_pt;
H = H0;
f = f_func(x);
g = g_func(x);
% d = length(init_pt);
% current_dis = 0;

% Line search
n = length(x);
fcn = @(n,x) deal(f_func(x),g_func(x));
% line search parameters
stp_init = 1;
ftol = 1e-2;
gtol = 5e-1;
% use the same setting as scipy 
% https://github.com/scipy/scipy/blob/3756215e9e51d3e5750cf2c7fe5bb014235b916a/scipy/optimize/_linesearch.py#L91
xtol = 1e-14;
stpmin = 1e-8;
stpmax = 1e3;
% stpmax = 1e4;
maxfev = 100;
for i_iter = 1:N_iter
    % More-Thuente line search for strong Wolfe condition
    p = -H*g;
    [x_new,f_new,g_new,stp,info,nfev] = cvsrch(fcn,n,x,f,g,p,stp_init,ftol,gtol,xtol, ...
                 stpmin,stpmax,maxfev);
    list_steps(i_iter) = nfev;
    %disp(stp);

    switch info
        case 1

        case 0
            disp('Improper input parameters')
        case 2
            disp('Relative width of the interval of uncertainty is at most xtol.')
        case 3
            disp('Number of calls to fcn has reached maxfev.')
        case 4
            disp('The step is at the lower bound stpmin.')
        case 5
            disp('The step is at the upper bound stpmax.')
        case 6
            disp('Rounding errors prevent further progress.')       
    end 
    % x_new = x-H*g;
    s = x_new-x;
    % g_new = g_func(x_new);
    y = g_new-g;
    rho = 1/(y'*s);
    % H = (eye(d)-rho*s*y')*H*(eye(d)-rho*y*s')+rho*(s*s');
    Hy = H*y;
    H = H-rho*(s*Hy'+Hy*s')+(rho^2*y'*Hy+rho)*(s*s');
    x = x_new;
    f = f_new;
    g = g_new;
    list_iter(:,1+i_iter) = x;
    if norm(s)<epsilon || norm(g)<epsilon
        break;
    end
end
list_iter = list_iter(:,1:1+i_iter);
list_steps = [0;list_steps(1:i_iter)];
end