Starting with\begin{align} q^{\prime }\left ( t\right ) +\frac{\left ( I_{1}-I_{3}\right ) }{I_{2}}\Omega \ r\left ( t\right ) & =0\tag{1}\\ r^{\prime }\left ( t\right ) +\frac{\left ( I_{2}-I_{1}\right ) }{I_{3}}\Omega \ q\left ( t\right ) & =0\tag{2} \end{align}
To decouple the ODE’s, take derivatives and substituting back, we find\begin{align} q^{\prime \prime }\left ( t\right ) +k\ q\left ( t\right ) & =0\tag{3}\\ r^{\prime \prime }\left ( t\right ) +k\ r\left ( t\right ) & =0\tag{4} \end{align}
Where \(k=\frac{\left ( I_{1}-I_{3}\right ) \left ( I_{2}-I_{2}\right ) }{I_{2}I_{3}}\Omega ^{2}\). The solution to the above is (for stability \(k>0\) )\begin{align} q\left ( t\right ) & =A\cos \sqrt{k}t+B\sin \sqrt{k}t\tag{5}\\ r\left ( t\right ) & =C\cos \sqrt{k}t+D\sin \sqrt{k}t\tag{6} \end{align}
Now, \(q\left ( 0\right ) =\frac{\text{Imp}}{I_{yy}}\), hence \(A=\frac{\text{Imp}}{I_{yy}}\) . To find \(B\) take derivative\begin{equation} q^{\prime }\left ( t\right ) =-\sqrt{k}\frac{\text{Imp}}{I_{yy}}\sin \sqrt{k}t+B\sqrt{k}\cos \sqrt{k}t\tag{7} \end{equation} But at \(t=0\) then we go back and use Eq. (1) to find \(q^{\prime }\left ( 0\right ) \), and equate the result to the above at \(t=0\). Eq (1), at \(t=0\), gives\[ q^{\prime }\left ( 0\right ) +\frac{\left ( I_{1}-I_{3}\right ) }{I_{2}}\Omega \ r\left ( 0\right ) =0 \] But \(r\left ( 0\right ) =0\,\), hence \(q^{\prime }\left ( 0\right ) =0\), and so this results in \(B=0\) in Eq.(7), hence the solution for \(q\left ( t\right ) \) is\[ q\left ( t\right ) =\frac{\text{Imp}}{I_{yy}}\cos \sqrt{k}t \] We do the same for \(r\left ( t\right ) .\) From Eq. (6) we find that \(C=0\) since \(r\left ( 0\right ) =0\), so now \(r\left ( t\right ) \) reduces to\[ r\left ( t\right ) =D\sin \sqrt{k}t \] Hence\begin{equation} r^{\prime }\left ( t\right ) =\sqrt{k}D\cos \sqrt{k}t\tag{8} \end{equation} To find \(D\), then we go back and use Eq. (2) to find \(r^{\prime }\left ( 0\right ) \), and equate the result to the above at \(t=0\). Eq (2), at \(t=0\), gives\[ r^{\prime }\left ( 0\right ) +\frac{\left ( I_{2}-I_{1}\right ) }{I_{3}}\Omega \ q\left ( 0\right ) =0 \] But \(q\left ( 0\right ) =\frac{\text{Imp}}{I_{yy}}\), hence\[ r^{\prime }\left ( 0\right ) =-\frac{\left ( I_{2}-I_{1}\right ) }{I_{3}}\Omega \ \frac{\text{Imp}}{I_{yy}}\] Therefore, equate the above to Eq. (8) evaluated at \(t=0\,\ \)gives\[ D=-\frac{\left ( I_{2}-I_{1}\right ) }{I_{3}\sqrt{k}}\Omega \ \frac{\text{Imp}}{I_{yy}}\] Hence Eq. (6) becomes\[ r\left ( t\right ) =\left ( -\frac{\left ( I_{2}-I_{1}\right ) }{I_{3}\sqrt{k}}\Omega \ \frac{\text{Imp}}{I_{yy}}\right ) \sin \sqrt{k}t \] so \(r\left ( t\right ) \) is sinusoidal as well.