% Discussed with: Gianmarco De Vita (MATLAB solver for determining \tau)

function pk = dogleg(B, g, deltak) 
    pnewton = - (B \ g);
    
    if norm(pnewton, 2) <= deltak
        pk = pnewton;
    else
        pu = - dot(g, g) / (g' * B * g) * g;
        
        if norm(pu, 2) > deltak
            pk = cauchy(B, g, deltak);
        else
            syms taux
            eqn = norm(pu + taux * (pnewton - pu))^2 == deltak^2;
            tau = solve(eqn, taux);
            tau = double(max(tau));
            pk = pu + tau * (pnewton - pu);
        end
    end
end