Eliminating $\vec N$ from these two equations we obtain \[ {d\vec L_G\over dt} = -\vec r_G \times M\left( {d^2 \vec r_G\over dt^2} +g\vec{e}_z \right) , \tag{2} \] which we must solve.
Expressing this in terms of the components of the angular velocity vector $\vec\omega$ in the body-fixed coordinate system ($XYZ$ system, where the $Z$ axis is aligned with the top’s axis), we obtain \[\left\{\begin{array}{rl}\displaystyle {d \omega_X\over dt} &\displaystyle = \left(1-{C\over A_0}\right)\omega_0\omega_Y + {Mr_Gg\over A_0}e_{zY} \\ \displaystyle {d \omega_Y\over dt} &\displaystyle = -\left(1-{C\over A_0}\right)\omega_0\omega_X - {Mr_Gg\over A_0}e_{zX} \\ \displaystyle {d \omega_Z\over dt} & = 0 \qquad\Rightarrow\quad \omega_Z = \omega_0\quad\mbox{(const.)} \end{array}\right. \tag{3} \] Here, $A$, $B (=A)$, and $C$ are the principal moments of inertia about the center of mass. $A_0 := A + M r_G^2$ is the moment of inertia about an axis through the contact point O, parallel to the $X$ or $Y$ axis. $\omega_0$ is the angular velocity around the axis of the top (i.e., $\omega_Z$), which is constant.
$e_{zX}$ and $e_{zY}$ are the $X$ and $Y$ components of the vertical unit vector $\vec e_z$ expressed in the body-fixed $XYZ$ coordinate system. Using the quaternion $R_q$ which transforms from the body frame to the lab frame, these may be obtained by \[ e_{zX}\,\hat i + e_{zY}\,\hat j + e_{zZ}\,\hat k = R_q^{-1}\hat k\, R_q . \tag{4} \] The time evolution of $R_q$ is given by \[ {d R_q\over dt} = {1\over 2}\ R_q\, \omega_q\,;\qquad \omega_q:= \omega_X \hat i + \omega_Y \hat j + \omega_Z \hat k . \tag{5} \] Solving the coupled system of equations (3)–(5) yields the dynamics the top.
In actual numerical computations, we eliminate the rotation in equation (3) by solving in a rotated coordinate frame ($X'Y'Z'$) that spins with angular velocity $-(1-C/A_0)\omega_0 \vec e_Z$ relative to the body frame ($XYZ$).
For details, see the computation notes(Jaanese).