We will provide a solution using a the Dyson expansion series. From the derivation we will see that the method can be easily extended to more complicated cases , ie when the power law function $\xi^{-2\beta}$ is replaced by a sum of different powers. For simplicity we only consider the homogenous ODE:
\begin{equation}
\frac{d^2 r_t}{d t^2} + U(t) r_t = 0
\end{equation}
where $U(t) = \omega^2/t^{2 \beta}$ .
That equation can be reduced to a following first order matrix ODE :
\begin{equation}
\frac{d \vec{x}_t}{d t} = \underbrace{\left( \begin{array}{rr} 0 && 1 \\ -U(t) && 0 \end{array} \right)}_{M(t)} \vec{x}_t
\end{equation}
where $\vec{x}_t := (r_t,v_t)$ where $v_t := \stackrel{.}{r}_t$. The solution to the ODE above, subject to initial conditions $r(t_0)=r_0$ and $v(t_0) = v_0$, is given by the Dyson series (time ordered exponential) as follows:
\begin{equation}
\vec{x}_t = \vec{x}_0 + \sum\limits_{p=1}^\infty
\underbrace{\int\limits_{t_0 \le \xi_{p-1} \le \cdots \le \xi_0 \le t} \left( \prod\limits_{l=0}^{p-1} M(\xi_l) d \xi_l \right)}_{G_p(t_0,t)} \cdot \vec{x}_0
\end{equation}
The terms in the product in the integrand are ordered in time, ie time increases from the left to the right. Now it is easy to check that the integrand reads:
\begin{eqnarray}
\prod\limits_{l=0}^{p-1} M(\xi_l) = (-1)^q \left\{
\begin{array}{cc}
\left(\begin{array}{cc} \prod\limits_{l=1}^q U(\xi_{2 l-1}) && 0 \\ 0 && \prod\limits_{l=1}^q U(\xi_{2l-2})\end{array}\right) && \mbox{if $p=2q$} \\
\left(\begin{array}{cc}0 && -\prod\limits_{l=1}^{q-1} U(\xi_{2l-1}) \\ \prod\limits_{l=1}^q U(\xi_{2l-2}) && 0\end{array}\right) && \mbox{if $p=2q-1$} \\
\end{array}
\right.
\end{eqnarray}
Now we use the result in Note 2 in Yet another multivariable integral over a simplex . and we write:
\begin{eqnarray}
G_p(t_0,t) = (-1)^q \left\{
\begin{array}{rr}
\left(\begin{array}{cc} \omega^q F_p^{\vec{\beta}_1}(t_0,t) && 0 \\ 0 && \omega^q F_p^{\vec{\beta}_2}(t_0,t)\end{array}\right)&& \mbox{if $p=2q$} \\
\left(\begin{array}{cc} 0 && -\omega^{q-1} F_p^{\vec{\beta}_3}(t_0,t)\\ \omega^q F_p^{\vec{\beta}_4}(t_0,t)&&0\end{array}\right)&& \mbox{if $p=2q-1$}
\end{array}
\right.
\end{eqnarray}
where
\begin{equation}
F_p^{\vec{\beta}}(t_0,t) := \sum\limits_{m=0}^p (-1)^m \frac{t_0^{-B_m+m} t^{-B_p+B_m+p-m}}{\prod\limits_{j=1}^p (j+B_{m-j}-B_m) \prod\limits_{j=1}^{p-m} (j+ B_m - B_{m+j})}
\end{equation}
where $B_j := \sum\limits_{l=1}^j \beta_j$ and the vectors $\vec{\beta}_l$ for $l=1,\cdots,4$ read:
\begin{eqnarray}
\vec{\beta}_1 &:= & \left(0,\beta,0,\beta,\cdots,0,\beta,0,\beta\right)_{p=2 q} \\
\vec{\beta}_2 &:= & \left(\beta,0,\beta,0,\cdots,\beta,0,\beta,0\right)_{p=2 q} \\
\vec{\beta}_3 &:= & \left(0,\beta,0,\beta,\cdots,0,\beta,0\right)_{p=2 q-1} \\
\vec{\beta}_4 &:= & \left(\beta,0,\beta,0,\cdots,\beta,0,\beta\right)_{p=2 q-1} \\
\end{eqnarray}
Therefore the final solution reads:
\begin{equation}
r_t = r_0 \sum\limits_{q=0}^\infty (-\omega)^q F_{2 q}^{\vec{\beta}_1}(t_0,t) - v_0 \sum\limits_{q=0}^\infty (-\omega)^q F_{2 q+1}^{\vec{\beta}_3} (t_0,t)
\end{equation}