Solution (source code)

= Solution

Solving the two momentum equations gives
$$
\dot\phi=\frac{L_z-L_3\cos\theta}{I_1\sin^2\theta},
\qquad
\dot\psi+\dot\phi\cos\theta=\frac{L_3}{I_3}.
$$
Substitution into the conserved energy reduces the <heavy symmetric top> to one-dimensional <effective potential> motion:
$$
\boxed{\frac12I_1\dot\theta^2+V_{\mathrm{eff}}(\theta)=E'},
$$
where
$$
\boxed{V_{\mathrm{eff}}(\theta)
=\frac{(L_z-L_3\cos\theta)^2}{2I_1\sin^2\theta}
+Mgl\cos\theta},
\qquad
\boxed{E'=E-\frac{L_3^2}{2I_3}}.
$$
This is the <heavy symmetric top reduction>.