function [xk, xs, gs] = BGFS(f, x0, H, max_itr, tol)
    syms x y;
    xk = x0;
    gf = gradient(f);
    
    fl = matlabFunction(f);
    gfl = matlabFunction(gf);
    
    tnorm = norm(at(gf, xk), 2);
    xs = zeros(size(x0, 1), max_itr + 1);
    xs(:, 1) = x0;
    gs = zeros(1, max_itr + 1);
    gs(1) = tnorm;
    k = 1;
    
    I = eye(2);
    
    while tnorm > tol && k <= max_itr 
        grad = gfl(xk(1), xk(2));
        p = -H * grad;
        
        alpha = backtracking(fl, gfl, p, xk, 1, 0.9, 1e-4);
        xk = xk + alpha * p;
        s = alpha * p;
        y = gfl(xk(1), xk(2)) - grad;
        
        rho = 1 / (y' * s);
        
        H = (I - (rho * s * y')) * H * (I - (rho * y * s')) ...
            + (1 / (y' * s) * (s * s'));

        k = k + 1;
        xs(:, k) = xk;
        tnorm = norm(gfl(xk(1), xk(2)), 2);
        gs(k) = tnorm;
    end
    xs = xs(:, 1:k);
    gs = gs(:, 1:k);
end