# Solution to a second order ordinary differential equation

Let $\beta >0$, $\gamma > 0$, $\omega > 0$ and $\xi >\xi_0$. The question is to show that the solution to the following inhomogenous ordinary differential equation: \begin{equation} \frac{1}{\omega^2} \frac{d^2 r_\xi}{d \xi^2} + \frac{1}{\xi^{2 \beta}} r_\xi = \xi^{-\gamma-1} \end{equation} subject to $r_{\xi_0} = r_0$ and $\dot{r}_{\xi_0}= \mu_0(-\omega)$ reads: \begin{eqnarray} (1) r_\xi =\\ C_\beta \sqrt{\xi} \left[ r_0(-\omega \xi_0^{1/2-\beta}) \left( J_{\frac{1}{2(-1+\beta)}}\left[x_\xi\right]J_{-1-\frac{1}{2(-1+\beta)}}\left[x_{\xi_0}\right]+ J_{\frac{-1}{2(-1+\beta)}}\left[x_\xi\right]J_{1+\frac{1}{2(-1+\beta)}}\left[x_{\xi_0}\right] \right) +\\ \mu_0(-\omega \xi_0^{1/2}) \left( -J_{\frac{1}{2(-1+\beta)}}\left[x_\xi\right]J_{-\frac{1}{2(-1+\beta)}}\left[x_{\xi_0}\right]+ J_{\frac{-1}{2(-1+\beta)}}\left[x_\xi\right]J_{\frac{1}{2(-1+\beta)}}\left[x_{\xi_0}\right] \right) +\\ \omega^2 \int\limits_{\xi_0}^\xi \eta^{-\gamma-1/2} \sum\limits_{s=\pm} s J_{\frac{s}{2(-1+\beta)}}\left[ x_\eta\right] J_{\frac{-s}{2(-1+\beta)}}\left[ x_\xi\right]d\eta \right] \end{eqnarray} where \begin{equation} C_\beta := \frac{\pi }{2(1-\beta) \sin\frac{\pi}{2(-1+\beta)}} \end{equation} and \begin{equation} x_\xi := \frac{\omega}{-1+\beta} \xi^{1-\beta} \end{equation}

where $J_\nu[x]$ is the Bessel function. Is it possible to simplify the integral over the Bessel functions?

Not sure if you'll find this useful or not, but if you set $\beta=1$ I found a special solution to your DEQ: $$r_\xi=\bigg(\frac{\omega^2}{\gamma^2-\gamma+\omega^2}\bigg)\xi^{1-\gamma}$$

• Thank you for that. If $\beta=1$ then one can show that the DEQ maps to an Euler equation and thus can be solved fairly easily. And as for that general solution, I have found it with help of Mathematica. Mathematica produced solutions for some specific values of beta and from that I deduced the general form. Then I checked it using the power series method. – Przemo Jun 27 '14 at 9:44
• Your special solution corresponds to the integral term in my generic solution. I have analyzed the $\beta \rightarrow 1$ limit of that term in the third answer from the top and I indeed retrieve your solution. – Przemo Jan 16 '15 at 18:51

We differentiate the function $r_\xi$ and we get: \begin{eqnarray} (2) \dot{r}_\xi = \\ C_\beta \omega \xi^{\frac{1}{2} -\beta}\left[ r_0(-\omega \xi_0^{1/2 - \beta})\left( J_{1+\frac{1}{2(-1+\beta)}}\left[x_\xi\right] J_{-1-\frac{1}{2(-1+\beta)}}\left[x_{\xi_0}\right] - J_{-1-\frac{1}{2(-1+\beta)}}\left[x_\xi\right] J_{1+\frac{1}{2(-1+\beta)}}\left[x_{\xi_0}\right] \right) +\\ \mu_0(-\omega \xi_0^{1/2})\left( -J_{1+\frac{1}{2(-1+\beta)}}\left[x_\xi\right] J_{-\frac{1}{2(-1+\beta)}}\left[x_{\xi_0}\right] - J_{-1-\frac{1}{2(-1+\beta)}}\left[x_\xi\right] J_{\frac{1}{2(-1+\beta)}}\left[x_{\xi_0}\right] \right)-\\ \omega^2 \int\limits_{\xi_0}^\xi \eta^{-\gamma-1/2} \sum\limits_{s=\pm} J_{\frac{s}{2(-1+\beta)}}\left[x_\eta\right] J_{-s-\frac{s}{2(-1+\beta)}}\left[x_\xi\right] \right] \end{eqnarray} Here I used the identity: \begin{equation} \frac{d}{d\xi} \sqrt{\xi} J_{\frac{s}{2(-1+\beta)}}\left[x_\xi\right]= s \omega \xi^{1/2-\beta} J_{s + \frac{s}{2(-1+\beta)}}\left[x_\xi\right] \end{equation}

Now we use equations (1) and (2) to check the initial conditions.

We set $\xi=\xi_0$ in equation one. It is readily seen that it is only the term proportional to $r_0$ that is non-zero. It reads: \begin{equation} (3) r_{\xi=\xi_0} = C_\beta \sqrt{\xi_0} r_0 (-\omega \xi_0^{1/2-\beta}) \cdot (-2) \frac{\sin(\pi \frac{-1}{2(-1+\beta)})}{\pi x_{\xi_0}} = r_0 \end{equation}

Now we set $\xi=\xi_0$ in equation (2).Likewise, it is readily seen that it is only the term proportional to $\mu_0$ that yields a non-zero contribution. It reads: \begin{equation} (4) \dot{r}_{\xi=\xi_0} = C_\beta \omega \xi_0^{1/2-\beta} \cdot \mu_0 (-\omega \xi_0^{1/2}) (-1) \frac{-2 \sin(\pi \frac{-1}{2(-1+\beta)})}{\pi x_{\xi_0}} = -\mu_0 \omega \end{equation} In equations (3) and (4) I used the identity: \begin{equation} J_{-a}\left[x\right] J_{1+a}\left[x\right] + J_{a}\left[x\right] J_{1-a}\left[x\right] = (-2) \frac{\sin(\pi a)}{\pi x} \end{equation}

Now that we see that the initial conditions are satisfied correctly. To finish the proof we need to differentiate equation (2) one more time. This requires a little bit more effort since now the integral term will produce an additional term. Yet the differentiations are pretty straightforward. I will complete this as soon as I find some more time.

We differentiate equation (2) one more time. We have: \begin{eqnarray} (3) \ddot{r}_\xi = \\ C_\beta \omega^2 \xi^{\frac{1}{2} -2\beta}\left[ r_0(-\omega \xi_0^{1/2 - \beta})\left( -J_{\frac{1}{2(-1+\beta)}}\left[x_\xi\right] J_{-1-\frac{1}{2(-1+\beta)}}\left[x_{\xi_0}\right] - J_{-\frac{1}{2(-1+\beta)}}\left[x_\xi\right] J_{1+\frac{1}{2(-1+\beta)}}\left[x_{\xi_0}\right] \right) +\\ \mu_0(-\omega \xi_0^{1/2})\left( J_{\frac{1}{2(-1+\beta)}}\left[x_\xi\right] J_{-\frac{1}{2(-1+\beta)}}\left[x_{\xi_0}\right] - J_{-\frac{1}{2(-1+\beta)}}\left[x_\xi\right] J_{\frac{1}{2(-1+\beta)}}\left[x_{\xi_0}\right] \right)-\\ \omega^2 \int\limits_{\xi_0}^\xi \eta^{-\gamma-1/2} \sum\limits_{s=\pm} s J_{\frac{s}{2(-1+\beta)}}\left[x_\eta\right] J_{-\frac{s}{2(-1+\beta)}}\left[x_\xi\right] \right]\\ - C_\beta \omega \xi^{1/2 - \beta}\cdot\omega^2 \xi^{-\gamma-1/2} \left( \sum\limits_{s=\pm} J_{\frac{s}{2(-1+\beta)}}\left[x_\eta\right] J_{-s-\frac{s}{2(-1+\beta)}}\left[x_\xi\right] \right) \end{eqnarray} Here we used the following identity: \begin{equation} \frac{d}{d\xi} \xi^{1/2 - \beta} J_{s + \frac{s}{2(-1+\beta)}}\left[x_\xi\right] = -s \omega \xi^{1/2 - 2 \beta} J_{\frac{s}{2(-1+\beta)}}\left[x_\xi\right] \end{equation} for $s=\pm$. By comparing equations (1) and (3) we see that the term in square brackets in (3) is just equal to $-\omega^2/\xi^{2\beta} r_\xi$ whereas the term outside the square bracket is just equal to $\omega^2 \xi^{-\gamma-1}$. Indeed, the term outside the square bracket reads: \begin{eqnarray} - C_\beta \omega \xi^{1/2 - \beta}\cdot\omega^2 \xi^{-\gamma-1/2} (-2) \frac{\sin\left(\pi \frac{-1}{2(-1+\beta)}\right)}{\pi x_\xi} =\\ -C_\beta \omega^2 \xi^{-\gamma-1} \cdot \frac{2(1-\beta)}{\pi} \sin\left( \pi \frac{-1}{2(-1+\beta)}\right)=\\ -C_\beta \omega^2 \xi^{-\gamma-1} \cdot \frac{-1}{C_\beta} = \omega^2 \xi^{-\gamma-1} \end{eqnarray}

Thus we have: \begin{equation} \ddot{r}_\xi = -\frac{\omega^2}{\xi^{2 \beta}} r_\xi + \omega^2 \xi^{-\gamma-1} \end{equation} which is equuivalent with the original ODE.

For simplicity we will provide an answer for the homogenous equation only. We also set $\omega=1$ In other words we seek to solve a following ODE: \begin{equation} \ddot{r}_t + \frac{1}{t^{2 \beta}} r_t =0 \end{equation} subject to $r_{t=t_0} = r_0$ and $\dot{r}_{t=t_0} = v_0$. By integrating this equation twice and solving the resulting integral equation iteratively it is easy to show that the general solution to the ODE reads: \begin{equation} r_t = r_0 J^{(t,t_0)} + v_0 I^{(t,t_0)} \end{equation} where \begin{eqnarray} J^{(t,t_0)} &:=& \sum\limits_{p=0}^\infty J_p^{(t,t_0)} \\ I^{(t,t_0)} &:=& -\sum\limits_{p=0}^\infty I_p^{(t,t_0)} \\ \end{eqnarray} and the quantities $J_p^{(t,t_0)}$ and $I_p^{(t,t_0)}$ are defined and calculated in closed form in Multivariable integral over a simplex . Now we calculate the fundamental solutions by doing the sum over $p$. We define $\theta := 1/(2\beta-2)$ and ${\mathcal C}_\beta := (\pi \theta)/ \sin(\pi \theta)$ and we have: \begin{eqnarray} &&-I{(t,t_0)}/ {\mathcal C}_\beta = \sum\limits_{s=0}^1\\ && \sum\limits_{p=0}^\infty \frac{(-1)^p}{(2 \beta-2)^{2 p} (p!)^2} \sum\limits_{q=0}^p \binom{p}{q} \binom{p}{q+(-1)^{1-s} \theta} t^{(2-2 \beta)(p-q)+s} t_0^{(2-2 \beta)q+1-s} (-1)^{1-s} =\\ && \sum\limits_{q_1=0}^\infty \sum\limits_{q_2=0}^\infty \frac{(-1)^{q_1} \left(\frac{t_0^{1-\beta}}{2(\beta-1)}\right)^{2 q_1}}{q_1!(q_1+(-1)^{1-s} \theta)!} \cdot \frac{(-1)^{q_2} \left(\frac{t^{1-\beta}}{2(\beta-1)}\right)^{2 q_2}}{q_2!(q_2+(-1)^{s} \theta)!} t^s t_0^{1-s} (-1)^{1-s} = \\ && \left(\frac{t_0^{1-\beta}}{2\beta-2}\right)^{(-1)^s \theta} \left(\frac{t ^{1-\beta}}{2\beta-2}\right)^{(-1)^{1-s} \theta} J_{(-1)^{1-s} \theta}\left(\frac{t_0^{1-\beta}}{\beta-1}\right) J_{(-1)^{ s} \theta}\left(\frac{t^{1-\beta}}{\beta-1}\right) t^s t_0^{1-s} (-1)^{1-s} = \\ && \sqrt{t t_0} (-1)^{1-s} J_{(-1)^{1-s} \theta}\left(\frac{t_0^{1-\beta}}{\beta-1}\right) J_{(-1)^{ s} \theta}\left(\frac{t^{1-\beta}}{\beta-1}\right) \end{eqnarray} In the second line from the top we used the closed form expressions for the quantities $I_p^{(t,t_0)}$. In the third line from the top we changed the summation variables to $(q_1,q_2)=(q,p-q)$. In the fourth line from the top we used the definition of Bessel functions and in the last line we simplified the result. Now, we do the same for the second fundamental solution: \begin{eqnarray} &&J{(t,t_0)}/ {\mathcal C}_\beta = \sum\limits_{s=0}^1\\ &&\sum\limits_{p=0}^\infty \frac{(-1)^p}{(2 \beta-2)^{2p-1} (p!)^2} \sum\limits_{q=0}^p \binom{p}{q} \binom{p}{q+(-1)^{1-s} \theta} t^{(2-2\beta)(p-q)+s} t^{(2-2\beta)q-s} (q - (1-s) \theta)(-1)^{1-s} = \\ &&(2 \beta-2) \sum\limits_{q_1=0}^\infty \sum\limits_{q_2=0}^\infty \frac{(-1)^{q_1} \left(\frac{t_0^{1-\beta}}{2(\beta-1)}\right)^{2 q_1}}{q_1!(q_1+(-1)^{1-s} \theta)!} \cdot \frac{(-1)^{q_2} \left(\frac{t^{1-\beta}}{2(\beta-1)}\right)^{2 q_2}}{q_2!(q_2+(-1)^{s} \theta)!} t^s t_0^{-s} (q_1-(1-s) \theta) (-1)^{1-s} \\ && (2\beta-2) \left(\frac{t^{1-\beta}}{2\beta-2}\right)^{2 s} \sum\limits_{q_1=0}^\infty \sum\limits_{q_2=0}^\infty \frac{(-1)^{q_1+s} \left(\frac{t_0^{1-\beta}}{2(\beta-1)}\right)^{2 q_1}}{q_1!(q_1+(-1)^{1-s} (\theta+1))!} \cdot \frac{(-1)^{q_2} \left(\frac{t^{1-\beta}}{2(\beta-1)}\right)^{2 q_2}}{q_2!(q_2+(-1)^{s} \theta)!} t^s t_0^{-s} (-1)^{1-s} = \\ && \frac{1}{t_0^\beta} \sqrt{t_0 t} (-1)^1 J_{(-1)^{1-s} (\theta+1)}\left(\frac{t_0^{1-\beta}}{\beta-1}\right) J_{(-1)^{ s} \theta}\left(\frac{t^{1-\beta}}{\beta-1}\right) \end{eqnarray} Again, in the second line from the top we used the closed form solution for the quantities $J_p^{(t,t_0)}$. In the third line from the top we changed the summation variables to $(q_1,q_2)=(q,p-q)$. In the last two remaining lines we firstly simplified the result and then expressed it through Bessel functions. The two fundamental solutions match to the homogenuous equation match those given in the formulation of the question.

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}

In here I am going to take the limit $\beta \rightarrow 1_+$ of the solution provided. For simplicity I analyze the last term, meaning the integral term, only. Using the asymptotic expansion for the Bessel functions: \begin{equation} \lim\limits_{|x| \rightarrow \infty}J_\theta (x) = \sqrt{\frac{2}{\pi x}} \cos(x - \frac{\pi}{2} \theta-\frac{\pi}{4}) \end{equation} we write the integral term as follows: \begin{eqnarray} &&\omega^2 C_\beta \sqrt{\xi} \cdot \left(\frac{\pi \omega}{2(\beta-1)}\right)^{-1}\cdot \\ &&\int\limits_{\xi_0}^{\xi} \eta^{-\gamma-1/2} \sum\limits_{s=\pm} s \cos(\omega\frac{ \eta^{1-\beta}}{\beta-1} - s \theta \frac{\pi}{2} - \frac{\pi}{4}) \cos(\omega\frac{ \xi^{1-\beta}}{\beta-1} + s \theta \frac{\pi}{2} - \frac{\pi}{4}) d\eta =\\ && \omega^2 \frac{\sqrt{\xi}}{\omega} \int\limits_{\xi_0}^{\xi} \eta^{-\gamma-1/2} \sin\left(\omega \frac{\eta^{1-\beta}-\xi^{1-\beta}}{\beta-1}\right) d\eta \\ &&\stackrel{=}{\beta \rightarrow 1_+} \omega \sqrt{\xi} \int\limits_{\xi_0}^{\xi} \eta^{-\gamma-1/2} \sin\left(\omega \log(\frac{\xi}{\eta})\right) d \eta = \\ && \frac{\omega^2}{ \omega^2 +(1/2- \gamma)^2}\left(\xi^{1-\gamma} -\frac{1}{\omega}\xi^{1/2} \xi_0^{1/2-\gamma} \sin(\omega \log(\frac{\xi}{\xi_0}) + \phi)\right) \end{eqnarray} where \begin{equation} \cos(\phi) = \frac{1/2-\gamma}{\sqrt{\omega^2 +(1/2- \gamma)^2}} \quad\mbox{and}\quad \sin(\phi) = \frac{\omega}{\sqrt{\omega^2 +(1/2- \gamma)^2}} \end{equation} Looking at this result we immediately recognize the special solution given in the first answer to this question.

Here we compute the integral over the Bessel functions. Note that this integral can be written as: \begin{equation} {\mathcal I} = \int\limits_{t_0}^t \left| \begin{array}{cc} f_+(\eta) & f_-(\eta) \\ f_+(\xi) & f_-(\xi) \end{array} \right| \frac{\omega^2 \eta^{-\gamma-1}}{{\mathcal W}[f_+,f_-](\eta)} d \eta \end{equation} where \begin{equation} f_\pm(\xi) := \sqrt{\xi} J_{\frac{\pm 1}{2 \beta-2}}\left[\omega \frac{\xi^{1-\beta}}{\beta-1}\right] \end{equation} are the two fundamental solutions to the homogeneous ODE and ${\mathcal W}[f_+,f_-](\xi)$ is the Wronskian of those solutions. From Abel's identity we see that the Wronskian does not depend on $\xi$. Therefore we can compute it at any value of $\xi$ in particular at $\xi \rightarrow \infty$. Using the asymptotic expression for the Bessel functions: \begin{equation} J_\alpha(z) = \sqrt{\frac{2}{\pi z}} \cos\left(z -\alpha \frac{\pi}{2} - \frac{\pi}{4}\right) \quad \mbox{for $z \rightarrow \infty$} \end{equation} a simple calculation shows that \begin{equation} {\mathcal W}[f_+,f_-](\xi) = \frac{\sin(\frac{\pi}{2 \beta-2})}{\frac{\pi}{2\beta -2}} = (C_\beta)^{-1} \end{equation} Now we write ${\mathcal I} = C_\beta \omega^2 {\mathcal J}_{-\gamma-1}$ where: \begin{eqnarray} &&{\mathcal J}_{-\gamma-1}= \left. \left| \begin{array}{cc} f_+(\eta) & f_-(\eta) \\ f_+(\xi) & f_-(\xi) \end{array} \right| \left(\frac{\eta^{-\gamma}}{-\gamma}\right) \right|_{\xi_0}^\xi + \int\limits_{\xi_0}^\xi \left| \begin{array}{cc} f_+^{'}(\eta) & f_-^{'}(\eta) \\ f_+(\xi) & f_-(\xi) \end{array} \right| \left(\frac{\eta^{-\gamma}}{\gamma}\right) d\eta \\ &&=\left. \left| \begin{array}{cc} f_+(\eta) & f_-(\eta) \\ f_+(\xi) & f_-(\xi) \end{array} \right| \left(\frac{\eta^{-\gamma}}{-\gamma}\right) \right|_{\xi_0}^\xi + \left. \left| \begin{array}{cc} f_+^{'}(\eta) & f_-^{'}(\eta) \\ f_+(\xi) & f_-(\xi) \end{array} \right| \left(\frac{\eta^{1-\gamma}}{(1-\gamma)\gamma}\right) \right|_{\xi_0}^\xi + \int\limits_{\xi_0}^\xi \left| \begin{array}{cc} f_+^{''}(\eta) & f_-^{''}(\eta) \\ f_+(\xi) & f_-(\xi) \end{array} \right| \left(\frac{\eta^{1-\gamma}}{(\gamma-1)\gamma}\right) d\eta \end{eqnarray} Here we have integrated by parts once and twice in the first and in the second line respectively. Now, since the functions $f_\pm$ satisfy the homogeneous ODE we get: \begin{eqnarray} &&{\mathcal J}_{-\gamma-1} + \frac{\omega^2}{(\gamma-1)\gamma} {\mathcal J}_{-\gamma+1-2 \beta} = \\ &&(C_\beta)^{-1} \frac{\xi^{1-\gamma}}{(\gamma-1) \gamma} + \left| \begin{array}{cc} f_+(\xi_0) & f_-(\xi_0) \\ f_+(\xi) & f_-(\xi) \end{array} \right| \frac{\xi_0^{-\gamma}}{\gamma} + \left| \begin{array}{cc} f_+^{'}(\xi_0) & f_-^{'}(\xi_0) \\ f_+(\xi) & f_-(\xi) \end{array} \right| \frac{\xi_0^{1-\gamma}}{(\gamma-1)\gamma} \end{eqnarray} The equation above is easily solved by recursion in $\gamma$. Assuming that $\beta > 1$ we have: \begin{eqnarray} &&{\mathcal I} = \omega^2 \frac{\xi^{1-\gamma}}{(\gamma)_{(2)}} F_{1,2} \left[ \begin{array}{rr} 1 \\ 1 + \frac{\gamma}{2\beta-2} & 1 + \frac{\gamma-1}{2\beta-2} \end{array}; -\frac{\omega^2}{4(\beta-1)^2 \xi^{2\beta-2}} \right] + \\ && C_\beta \omega^2 \frac{\xi_0^{-\gamma}}{\gamma} F_{1,2} \left[ \begin{array}{rr} 1 \\ 1 + \frac{\gamma}{2\beta-2} & \frac{\gamma-1}{2\beta-2} \end{array}; -\frac{\omega^2}{4(\beta-1)^2 \xi_0^{2\beta-2}} \right]\left(f_+(\xi_0) f_-(\xi) - f_-(\xi_0) f_+(\xi)\right) + \\ && C_\beta \omega^2 \frac{\xi_0^{1-\gamma}}{(\gamma)^{(2)}} F_{1,2} \left[ \begin{array}{rr} 1 \\ 1 + \frac{\gamma}{2\beta-2} & 1+\frac{\gamma-1}{2\beta-2} \end{array}; -\frac{\omega^2}{4(\beta-1)^2 \xi_0^{2\beta-2}} \right]\left(f_+^{'}(\xi_0) f_-(\xi) - f_-^{'}(\xi_0) f_+(\xi)\right) \end{eqnarray}