Solution (source code)

= Solution

Put $L_z=L_3(1-\epsilon^2/2)$ and seek $\theta=O(\epsilon)$. Using $\cos\theta=1-\theta^2/2+O(\theta^4)$ and $\sin^2\theta=\theta^2+O(\theta^4)$, the $\theta$-dependent part of the effective potential is
$$
V_{\mathrm{eff}}(\theta)
=\text{constant}
+\frac{L_3^2-4MglI_1}{8I_1}\theta^2
+\frac{L_3^2\epsilon^4}{8I_1\theta^2}
+O(\epsilon^4).
$$
It has a local minimum precisely when
$$
\boxed{L_3^2>4MglI_1},
$$
which is the <gyroscopic stabilization of an inverted symmetric top> condition. Differentiating the displayed approximation gives
$$
\theta_{\mathrm{eq}}^4
=\frac{L_3^2\epsilon^4}{L_3^2-4MglI_1},
$$
and hence
$$
\boxed{\theta_{\mathrm{eq}}
=\epsilon\left(\frac{L_3^2}{L_3^2-4MglI_1}\right)^{1/4}}.
$$
Taking $L_3>0$ without loss of orientation, the steady <precession> rate is
$$
\dot\phi
=\frac{L_z-L_3\cos\theta_{\mathrm{eq}}}
{I_1\sin^2\theta_{\mathrm{eq}}}
=\frac{L_3}{2I_1}
\left(1-\frac{\epsilon^2}{\theta_{\mathrm{eq}}^2}\right)+O(\epsilon^2),
$$
so
$$
\boxed{\dot\phi
=\frac{L_3-\sqrt{L_3^2-4MglI_1}}{2I_1}
+O(\epsilon^2)}.
$$