My colleague and I are trying to understand how to solve the following set of ODEs with non-constant coefficients:
$$ c_1\left[f''(x)+\frac{1}{x}f'(x)-\left(a_1^2+\frac{1}{x^2}\right)f(x)\right]+c_2\left[ g''(x)+\frac{1}{x}g'(x)-\left(a_2^2+\frac{1}{x^2}\right)g(x)\right]=0\,, $$
$$ c_3\left[f''(x)+\frac{1}{x}f'(x)-\left(a_3^2+\frac{1}{x^2}\right)f(x) \right]+c_4\left[ g''(x)+\frac{1}{x}g'(x)-\left(a_4^2+\frac{1}{x^2}\right)g(x)\right]=0\,. $$
Here, $c_i\in\mathbb{R}$ and $a_i^2\in\mathbb{R}$, $i=1,\cdots,4$ are real-valued constants. In this paper, the following Ansatz is made to solve the problem
$$ f(x)=h_1 I_1(\lambda_1 x)+h_2 I_1(\lambda_2 x), \qquad g(x)=k_1 I_1(\lambda_1 x)+k_2 I_1(\lambda_2 x), $$
where $I_1$ denotes the modified Bessel function of the first kind and of order one and $h_i$ and $k_i$ are constants. Moreover, $\lambda_{1,2}$ are the positive solutions of
$$ (\lambda_{1,2}^2-a_1^2)(\lambda_{1,2}^2-a_1^2)c_1c_4-(\lambda_{1,2}^2-a_2^2)(\lambda_{1,2}^2-a_3^2)c_2c_3=0. $$
We understand the steps of the authors of the paper to some extent, but would like to know how we can be sure that the solution is unique. We think that $f$ and $g$ should have two linear independent solutions because the ODEs are second order. It also makes sense to use Bessel functions in the Ansatz because each bracket in the first two equations can be transformed to a modified Bessel differential equation. Since regularity for $x=0$ should be guaranteed, modified Bessel functions of the second kind are dropped. However, why should the coordinate transformations of the argument of the Bessel functions used for $f$ and $g$ be the same? Couldn't, for example, $g(x)=k_1 I_1(\lambda_3 x)+k_2 I_1(\lambda_4 x)$ be also used?
We would be grateful for some comments.