Rolling Disk Constraint
The disk is described by lateral position $(x,y)$ and Euler angles $(\theta,\phi,\psi)$ (the dynamic variables), and constants of radius $r$ and mass per unit area $\lambda$. The Lagrangian (calculated through Mathematica) is:
\begin{aligned} L&=\frac{1}{16} \pi \lambda r^2 \Big(-16 g r \sin (\phi )+r^2 \Big(8 \dot{\theta} \dot{\psi} \cos (\phi )\\ &+\dot{\theta} ^2 (\cos (2 \phi )+3)+4 \dot{\psi} ^2+\dot{\phi} ^2 (4 \cos (2 \phi )+6)\Big)\\ &+8 \dot{x}^2+8 \dot{y}^2\Big) \end{aligned}With constraints:
\begin{array}{l} r \sin (\theta ) d\psi +r d\theta \sin (\theta ) \cos (\phi )+r \cos (\theta ) d\phi \sin (\phi )+dx =0\\ -r \cos (\theta ) d\psi +r \sin (\theta ) d\phi \sin (\phi )-r d\theta \cos (\theta ) \cos (\phi )+dy =0\\ \end{array}See my stackexchange thread here for more detail.
After we apply the method of undetermined multipliers to eliminate dependence on x and y, we get the final differential equations:
\begin{array}{l} \ddot{\theta}=2 \dot{\psi} \dot{\phi} \csc (\phi ) \\ \ddot{\phi}=-\frac{4 g \cos (\phi )}{5 r}-\frac{1}{5} \dot{\theta} \sin (\phi ) (5 \dot{\theta} \cos (\phi )+6 \dot{\psi} ) \\ \ddot{\psi}= \dot{\phi} \left(\frac{5}{3} \dot{\theta} \sin (\phi )-2 \dot{\psi} \cot (\phi )\right) \\ \end{array}Here you can see the numerical solution animated: