Solution (source code)

= Solution

Eliminate the cyclic velocities using $\dot\phi=(L_z-L_3\cos\theta)/(I_1\sin^2\theta)$. The conserved <energy> becomes
$$
E-\frac{L_3^2}{2I_3}=\frac{I_1}2\dot\theta^2+V_{\rm eff}(\theta),\qquad\boxed{V_{\rm eff}=\frac{(L_z-L_3\cos\theta)^2}{2I_1\sin^2\theta}+Mgl\cos\theta.}
$$
Differentiating the <effective potential> gives
$$
-\frac{dV_{\rm eff}}{d\theta}=I_1\dot\phi^2\sin\theta\cos\theta-L_3\dot\phi\sin\theta+Mgl\sin\theta,
$$
which is also the $\theta$ <Euler-Lagrange equation>. Hence \b[$I_1\ddot\theta=-V_{\rm eff}'(\theta)$]. Deriving it directly from the <Euler-Lagrange equation> ensures it also holds at turning points, where division by $\dot\theta$ in an <energy> argument would be invalid.