I have obtained the following coupled two ODEs as the subsidiary equations of two coupled PDEs in Laplace domain. I was trying to solve them by decoupling them and getting a higher order differential equation form for each ODE but the solution later became so complicated. I wonder if there is a method to solve them.
First ODE. $$ y_1''(x) + \frac{C_1}{x}y_1'(x)= C_2 y_1(x) + C_3 y_2(x)$$
Second ODE. $$ y_2''(x) + \frac{C_4}{x}y_2'(x)= C_5 y_2(x) + C_6 y_1(x)$$
For: First ODE: $$ y_1''(x) + \frac{C_1}{x}y_1'(x)= C_2 y_1(x) + C_3 y_2(x)$$ Second ODE: $$ y_2''(x) + \frac{C_4}{x}y_2'(x)= C_5 y_2(x) + C_6 y_1(x)$$
proceed as follows: multiply the first by $C_{6}$ and the second by $C_{2}$ and subtract to obtain $$C_{6} \, \left(y_{1}^{''} + \frac{C_{1}}{x} \, y_{1}^{'} \right) = \lambda = C_{2} \left( y_{2}^{''} + \frac{C_{4}}{x} \, y_{2}^{'} \right) + (C_{6} C_{3} - C_{2} C_{5}) y_{2},$$ where $\lambda$ is a separation constant. Now, \begin{align} y_{1}^{''} + \frac{C_{1}}{x} \, y_{1}^{'} &= \frac{\lambda}{C_{6}} \\ y_{2}^{''} + \frac{C_{1}}{x} \, y_{2}^{'} + \frac{C_{6}C_{3} - C_{2}C_{5}}{C_{2}} \, y_{2} &= \frac{\lambda}{C_{2}}. \end{align}
The $y_{1}$ equation can be solved as follows: \begin{align} y_{1}^{''} + \frac{C_{1}}{x} \, y_{1}^{'} &= \frac{\lambda}{C_{6}} \\ \frac{1}{x^{C_{1}}} \, \frac{d}{dx} \left( x^{C_{1}} \, y_{1}^{'} \right) &= \frac{\lambda}{C_{6}} \\ x^{C_{1}} \, \frac{d y_{1}}{dx} &= \frac{\lambda \, x^{C_{1} + 1}}{C_{6} \, (C_{1} + 1)} + d_{0} \\ \frac{d y_{1}}{dx} &= \frac{\lambda \, x}{C_{6} \, (C_{1} + 1)} + \frac{d_{0}}{x^{C_{1}}}\\ y_{1} &= \frac{\lambda \, x^{2}}{2 \, C_{6} \, (C_{1} + 1)} - \frac{d_{0}}{(C_{1} -1) \, x^{C_{1}-1}} + d_{1}. \end{align}
The second equation is solvable in terms of Bessel functions as seen by
$$f'' + \frac{a}{x} \, f' + b f = c$$
has the solution
$$f(x) = x^{(1-a)/2} \, \left(A_{0} \, J_{\beta}(\sqrt{b} \, x) + B_{0} \, Y_{\beta}(\sqrt{b} \, x) \right) + \frac{c}{b},$$
where $\beta = (1-a)/2$.
Given the conditions for the equations and the solution then the forms of solutions for $y_{1}$ and $y_{2}$ could be reduced.
What's your background? I'm not sure how much detail to include.
================================
Method One Note you can isolate y2(x) in the first ODE. Take the corresponding derivatives and plug them into the second ODE and you get a new ODE completely in terms of y1(x).
==========================
Method 2 Use the Method of Frobenius where y1 is expressed in terms of one set of coefficients and y2 is expressed as a series with a different set of coefficients, then proceed as normal. http://mathworld.wolfram.com/FrobeniusMethod.html
=============================================
Method 3, I'm a bit fuzzy on this one, but in p Solve the equation you get if you assume C3=0 and substitute y0 in place of y1: y0''+C1/x * y0' = C2y0. Again you should be able to get this into Bessel's Equation. Let y2=v*y1, swapping into both equations. Simplify. Then swap y1=u*y0. You should get simpler equations.
2 potential ideas:
(1) If you're interested in the behavior away from $x=0$:
Using the usual reduction of order substitutions $y_1'=y_3$ and $y_2'=y_4$, we can rewrite as a first order linear system $Y'(x) = A(x)Y(x)$, where
$$ Y(x) = \begin{pmatrix}y_1(x)\\y_2(x)\\y_3(x)\\y_4(x)\end{pmatrix}, \ Y'(x) = \begin{pmatrix}y_1'(x)\\y_2'(x)\\y_3'(x)\\y_4'(x)\end{pmatrix}, \text{ and } A(x) = \begin{pmatrix}0&0&1&0 \\ 0&0&0&1 \\ C_2&C_3&-C_1x^{-1}&0 \\ C_6&C_5&0&-C_4x^{-1} \end{pmatrix} $$
Let $B(x) = \int_{x_0}^x A(s)\text ds$ for the initial point $x_0$. Then, if $AB = BA$ for all $x$, $Y(x) = \exp(B(X))$ is the solution.
(2) If you're interested in the behavior near $x=0$:
Let $M(x) = xA(x)$, and rewrite the equation as $xY'(x) = M(x)Y(x)$. Now, $M(x)$ is holomorphic at $x=0$ and we can proceed with the method laid out in Chapter II of Wasow's "Asymptotic expansions for ordinary differential equations" $-$ essentially a matrix version of Frobenius/power series method.
