function [xk, xs, gnorms] = trust_region(f, delta_hat, delta0, eta, ...
    x0, tol, max_n)
    xs = zeros(2, max_n);
    gnorms = zeros(max_n);
    
    xk = x0;
    deltak = delta0;
    
    fl = vecLambda(matlabFunction(f));
    gl = vecLambda(matlabFunction(gradient(f)));
    hl = vecLambda(matlabFunction(hessian(f)));
    
    xs(:, 1) = x0;

    for i = 2:max_n
        fk = fl(xk);
        B = hl(xk);
        g = gl(xk);
        
        gnorms(i - 1) = norm(g);
       
        if gnorms(i - 1) < tol
            gnorms = gnorms(1:i-1);
            xs = xs(:, 1:i-1);
            break
        end

        pk = dogleg(B, g, deltak);
        
        rho_k = (fl(xk) - fl(xk + pk)) / ...
            (qf(B, g, [0;0], fk) - qf(B, g, pk, fk));
        
        if rho_k < 1/4
            deltak = 1/4 * deltak;
            
        % When comparing the method's execution with some classmates, we
        % found some numerical instability in the comparison between the
        % step's norm and the trust region radius. To sum up, using
        % different Matlab versions (2020b vs 2019a) we obtained different
        % results on that comparison (true vs. false) on seemingly
        % identical values. It seems there is some difference w.r.t.
        % comparisons in the unnormalized double range, and therefore we
        % both approximated that equality with the subtraction you will find
        % below
        elseif rho_k > 3/4 && (norm(pk, 2) - deltak) < eps
            deltak = min(2 * deltak, delta_hat);
        end
        % otherwhise do not change delta
        
        if rho_k > eta
            xk = xk + pk;
        end
        % otherwise do not change xk   
        
        xs(:, i) = xk;
    end
end

% Convert lambda to accept vector parameters
function vl = vecLambda(fl) 
    vl = @(x) fl(x(1), x(2));
end

% Compute quadratic form
function y = qf(B, g, p, fk)
    y = fk + 1/2 * p' * B * p + dot(g, p);
end