Bessel ODE 1
The ordinary differential equation by Bessel is defined by:
$$
x^2\, y_\lambda'' + x\, y_\lambda' + (x^2 - \lambda^2)\, y_\lambda = 0
$$
Everything in the real numbers. With operator calculus it reads:
$$
\left[x^2\frac{d^2}{dx^2} + x \frac{d}{dx} + x^2-\lambda^2 \right]\, y_\lambda(x) = 0
$$
This can be simplified using:
$$
x^2\frac{d^2}{dx^2} + x \frac{d}{dx} = \left(x\frac{d}{dx}\right)^2
\quad \Longrightarrow \quad
\left[\left(x\frac{d}{dx}\right)^2 + x^2-\lambda^2 \right]\, y_\lambda(x) = 0
$$
Let's try to decompose the operator in factors:
$$
\left(x\frac{d}{dx}\right)^2 + x^2-\lambda^2 =
\left[x\frac{d}{dx} + (ax + b)\right]\left[x\frac{d}{dx} - (ax + b)\right] = \\
\left(x\frac{d}{dx}\right)^2 - (a^2x^2+2abx+b^2) - ax
$$
It follows that:
$$
-a^2 = 1 \quad \Longrightarrow \quad a = \pm\, i \\
2ab+a = 0 \quad \Longrightarrow \quad b = -\frac{1}{2} \\
b^2 = \lambda^2 \quad \Longrightarrow \quad \lambda = \pm \frac{1}{2}
$$
We've got two values for $\lambda$. However, it's trivial that if
$y_\lambda(x)$ is a solution, then $y_{-\lambda}(x)$ is a solution as well
(bacause $\lambda$ does only occur as its square $\lambda^2$).
So, without loss of generality, we can restrict ourselves to solutions
with $\lambda \ge 0$.
The differential equation may now be written as follows:
$$
\left[x\frac{d}{dx} + i\,x - \frac{1}{2}\right]
\left[x\frac{d}{dx} - i\,x + \frac{1}{2}\right]\, y_\lambda(x) = 0\\
x\left[\frac{d}{dx} + i - \frac{1}{2x}\right]
x\left[\frac{d}{dx} - i + \frac{1}{2x}\right]\, y_\lambda(x) = 0
$$
With the well known formula in our non-volatile memory
- it's at the bottom line of the general theory:
$$
\frac{d}{dx}\pm\left(i-\frac{1}{2x}\right) =
e^{\mp\int\left(i-\frac{1}{2x}\right)dx} \frac{d}{dx}
e^{\pm\int\left(i-\frac{1}{2x}\right)dx} =\\
e^{\mp ix}e^{\pm\frac{1}{2}\ln(x)}\frac{d}{dx}e^{\pm ix}e^{\mp\frac{1}{2}\ln(x)}=
e^{\mp ix}x^{\pm\frac{1}{2}}\frac{d}{dx}e^{\pm ix}x^{\mp\frac{1}{2}}
$$
Giving for the differential equation:
$$
x\,e^{-ix}x^{\frac{1}{2}}\frac{d}{dx}e^{ix}x^{-\frac{1}{2}}
x\,e^{ix}x^{-\frac{1}{2}}\frac{d}{dx}e^{-ix}x^{\frac{1}{2}}\, y_\lambda(x) = 0
\quad \Longrightarrow \\
e^{2ix}\frac{d}{dx}e^{-ix}x^{\frac{1}{2}}\, y_\lambda(x) = c_1
\quad \Longrightarrow \quad
e^{-ix}x^{\frac{1}{2}}\, y_\lambda(x) = c_1 e^{-2ix} + c_2
$$
Leading to the solution, for $\lambda = \frac{1}{2}$ :
$$
y_{\frac{1}{2}}(x) = \frac{c_2 e^{+ix} + c_1 e^{-ix}}{\sqrt{x}}
= \frac{C_1 \cos(x) + C_2 \sin(x)}{\sqrt{x}}
$$
With $C_1,C_2$ real valued, because we have said: Everything in the real numbers.