On the right hand side you have $\theta$, $\dot\theta$, $r$ and $\dot r$. They should be $\theta_1$, $\theta_2$, $r_1$ and $r_2$.
Differentiate the first wrt $t$ to gain an expression for $\ddot y$:
$\dddot{x} = \frac{\omega_1^2}{2} \dot x + \omega_2 \ddot{y}$
Substitute $\ddot{y} = \frac{\omega_1^2}{2} y - \omega_2 \dot{x}$ to get:
$\dddot{x} = \frac{\omega_1^2}{2} \dot x + \frac{\omega_1^2 \omega_2}{2} y - \omega_2^2 \dot{x}$
Rearrange: $\frac{\omega_1^2 \omega_2}{2} y =\dddot{x} - \frac{\omega_1^2}{2} \dot x + \omega_2^2 \dot{x}$
Differentiate: $\frac{\omega_1^2 \omega_2}{2} \dot y =\ddddot{x} + \frac{2\omega_2^2 -\omega_1^2}{2} \ddot x$
Recall that $\ddot{x} = \frac{\omega_1^2}{2} x + \omega_2 \dot{y} \Rightarrow \omega_2 \dot{y}=\ddot{x} - \frac{\omega_1^2}{2} x$
Thus: $\frac{\omega_1^2}{2} \left (\ddot{x} - \frac{\omega_1^2}{2} x \right ) =\ddddot{x} + \frac{2\omega_2^2 -\omega_1^2}{2} \ddot x$
... which becomes $\ddddot{x} + \frac{2\omega_2^2 -\omega_1^2}{2} \ddot x - \frac{\omega_1^2}{2} \left (\ddot{x} - \frac{\omega_1^2}{2} x \right )=0 $
or $\ddddot{x} + \left(\omega_2^2 -\omega_1^2 \right ) \ddot x - \frac{\omega_1^4}{2}x =0 $
Auxiliary equation $\lambda^4+p\lambda^2-q=0$ where $p=\left(\omega_2^2 -\omega_1^2 \right )$ and $q=\omega_1^4$
$\lambda^2={-p+\sqrt{p^2+4q} \over 2}$ or $\lambda^2={-p-\sqrt{p^2+4q} \over 2}$
$\lambda_1=\sqrt{{-p+\sqrt{p^2+4q} \over 2}}$
$\lambda_2=-\sqrt{{-p+\sqrt{p^2+4q} \over 2}}$
$\lambda_3=\sqrt{{-p-\sqrt{p^2+4q} \over 2}}$
$\lambda_4=-\sqrt{{-p-\sqrt{p^2+4q} \over 2}}$
Best Answer
Let:
$z = x_1$
$z' = x_1' = x_2$
$z'' = x_1'' = x_2' = -\dfrac{-GM_E}{(r_0+x_1)^2} + (r_0+x_1)(\omega_0+x_4)^2 + \dfrac{F_r}{m}$
$\phi = x_3$
$\phi' = x_3' = x_4$
$\phi'' = x_3'' = x_4' = -2 \dfrac{(\omega_0+x_4)x_2}{r_0+x_1}+\dfrac{F_\theta}{m(r_0+x_1)}$
So now, we have a system of four-first order ordinary differential equations with $x_1', x_2', x_3', x_4'$.