function tau_max = tauMaxSolver(tau, L, Lmax, beta, m)

    % RHS dell'invariante
    R = (beta*tau)^m * ( L + (m/(m+1))*beta*tau );

    % funzione g(z) che sarebbe g(Tau_max) con Tau_max che è la nostra
    % incognita
    g = @(z) (beta*z)^m * ( Lmax + (m/(m+1))*beta*z ) - R;

    % bounds iniziali
    z_low = 0;              %limite fisico di Tau_max e Tau
    z_high = max(tau,1);    % così impongo già che Tau_max sia maggiore di Tau

    % aumento z_high finché g(z_high)>0 per garantire l'ipotesi del teorema
    % degli zeri
    while g(z_high) < 0
        z_high = 2*z_high;
    end

    % solve implicito con fzero che è una function di matlab e banalmente
    % risolve l'equazione trovando lo zero quindi trovando il valore di
    % Tau_max
    tau_max = fzero(g,[z_low z_high]);

end