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'$.
You have a coupled system of first order ODEs
\begin{align}
x' &= p \quad (1) \\
p' &= -x \quad (2)
\end{align}
This can be turned into two second order ODEs by differentiating both $x'$ and $p'$ again i.e
\begin{align}
x'' &= p' \quad \text{from (1)} \\
&= -x \quad \text{from (2)} \\\\
p'' &= -x' \quad \text{from (2)} \\
&= -p \quad \text{from (1)}
\end{align}
The boundary conditions can also be derived using your ODE system and the initial conditions for the PDE problem
\begin{align}
x(s,0) &= s, \quad x'(s,0) = 1 \\
p(s,0) &= 1, \quad p'(s,0) = -s
\end{align}
They are the same ODE, so we will just solve the one in $x$ here. There are two ways to solve this apart from using the matrix formulation. First, the ODE is linear so making an ansatz $x = e^{\lambda t}$ gives us $\lambda^{2} = -1 \implies \lambda = \pm i$, hence
\begin{align}
x &= Ae^{it} + Be^{-it} \\
&= A(\cos(t) + i\sin(t)) + B(\cos(t) - i \sin(t)) \\
&= C_{1} \cos(t) + C_{2} \sin(t)
\end{align}
where $C_{1} = A + B, C_{2} = i(A - B)$. The second way is to multiply the ODE by $x'$ and integrate i.e
\begin{align}
x' x'' &= -x x' \\
\implies \frac{1}{2} x'^{2} &= - \frac{1}{2} x^{2} + k_{1} \\
\implies x' &= \pm \sqrt{k_{1} - x^{2}} \\
\implies \int \frac{dx}{\sqrt{k_{1} - x^{2}}} &= \pm \int dt \\
\implies \arcsin \left( \frac{x}{k_{1}} \right) &= \pm t + k_{2} \\
\implies x &= k_{1} \sin(k_{2} \pm t) \\
&= c_{1} \cos(t) + c_{2} \sin(t)
\end{align}
where $c_{1} = k_{1} \sin(\pm k_{2}), c_{2} = k_{1} \cos(\pm k_{2})$. Now, applying your conditions we find
\begin{align}
x(0) &= s \\
&= C_{1} \cos(0) + C_{2} \sin(0) \\
&= C_{1} \\
x'(0) &= 1 \\
&= -s \sin(0) + C_{2} \cos(0) \\
&= C_{2}
\end{align}
Hence the particular solution is given by $x(s,t) = s \cos(t) + \sin(t)$. The same approach is used to solve the ODE in $p$.
Best Answer
This line is not correct you inverted $u,v$: $$\left(\begin{matrix}\dot u(t)\\\dot v(t)\\ \end{matrix} \right)=\left(\begin{matrix}0&1\\\frac{1}{1-t}&\frac{t}{t-1}\\ \end{matrix} \right)\left(\begin{matrix}v(t)\\u(t)\\ \end{matrix} \right)$$ It should be: $$\left(\begin{matrix}\dot u(t)\\\dot v(t)\\ \end{matrix} \right)=\left(\begin{matrix}0&1\\\frac{1}{1-t}&\frac{t}{t-1}\\ \end{matrix} \right)\left(\begin{matrix}u(t)\\v(t)\\ \end{matrix} \right)$$ So that the second order DE is: $$x''(t)-\frac{t}{t-1}x'(t)-\frac{1}{1-t}x(t)=0$$
Edit
@psyph Applying Liouville's formula I got this general solution:
$$|\phi (t)|=c_1 \exp \int \frac {t}{t-1} dt$$ $$|\phi (t)|=c_1e^t(t-1)$$ We also have that: $$|\phi (t)|=u(t)e^t-v(t)e^t=e^t(u(t)-v(t))$$ $$\implies u(t)e^t-v(t)e^t=c_1e^t(t-1)$$ $$ u(t)-v(t)=c_1(t-1)$$ from the original equation we have: $$u'(t)=v(t) \implies u'(t)-u(t)=-c_1(t-1)$$ using method of integrating factor $$ (u(t)e^{-t})'=-c_1(t-1)e^{-t}$$ Integrate: $$ u(t)e^{-t}=-c_1(-(t-1)e^{-t}-e^{-t})+c_2$$ $$ u(t)=-c_1(-(t-1)-1)+c_2e^{t}$$ $$ \implies u(t)=tc_1+c_2e^{t}$$ And for $v(t)$: $$v(t)=u'(t)=c_1+c_2e^{t}$$ Finally : $$\left(\begin{matrix} u(t)\\\ v(t)\\ \end{matrix} \right)=\left(\begin{matrix}c_1t+c_2e^t\\c_1+c_2e^t\\ \end{matrix} \right)$$ $$\left(\begin{matrix} u(t)\\\ v(t)\\ \end{matrix} \right)=c_1\left(\begin{matrix}t\\1\\ \end{matrix} \right)+c_2\left(\begin{matrix}e^t\\e^t\\ \end{matrix} \right)$$