Approximate analytic solution of Schroedinger equation with arbitrary power potential I'm solving the following Schroedinger equation in the domain $r>0$
$\psi''(r) + \left(E-\frac{a}{r^b}\right)\psi(r)=0 $,
where $0 < b < 2$ and $a, E$ are positive constants. Primarily I'm interested in the asymptotical power behavior of the solution as $r\to 0$. To be complete in the description of the problem, I fix my boundary conditions at $r\to +\infty$ as a plane wave ansatz.
I did a lot of DSolve with Mathematica and found out that $\psi(r)\to const$ as $r\to 0$. It gave me a hint for the power series solution, that one of the terms in the expansion
$\psi(r) = \sum\limits_i a_i r^{\alpha_i}$ 
for some (noninteger) $\alpha_i>0$ should cancel with $a/r^b$. However, even with this assumption there are a bunch of terms which do not cancel. Power expansion does not seem to work in this case as well as the WKB approximation (double checked numerically). What are the other known methods to find an approximate asymptotic behavior in this case? A good reference would be great too!
 A: As $r$ approaches zero, the coefficient of $\psi$ becomes dominated by the contribution of $-a/r^b$. This means that in the vicinity of zero your solution is dominated by solution of the following equation
$$
\psi''-\frac{a}{r^b}\psi=0.
$$
This can be demonstrated more rigorously (and, also, refined to the higher accuracy) by scaling into the vicinity of $r=0$. This equation possesses an explicit solution in terms of the modified Bessel functions:
$$
\psi=C_1\sqrt{r}K_{\frac{1}{2-b}}\left(\frac{2\sqrt{a}}{b-2}r^{1-b/2}\right)+C_2\sqrt{r}I_{\frac{1}{2-b}}\left(\frac{2\sqrt{a}}{b-2}r^{1-b/2}\right).
$$
The latest edition of NIST Handbook tells us that for $z\to 0$
$$
I_{\nu}(z)\sim (\frac{1}{2}z)^{\nu}/\Gamma(\nu+1)
\quad\mbox{and}\quad
K_{\nu}(z)\sim \frac{1}{2}\Gamma(\nu)(\frac{1}{2}z)^{-\nu}
$$
so the actual constant at $r\to0$ depends only on $C_1$, unless $C_1$ vanishes. It seems pretty clear that there are no singularities here; typically you should expect that
$$
\psi\sim \frac{C_1}{\Gamma(\frac{3-b}{2-b})}\left(\frac{\sqrt{a}}{b-2}\right)^{\frac{1}{2-b}}
\quad\mbox{for}\quad r\to 0.
$$
Specific values of $C_1$ and $C_2$ are less trivial to obtain, but they can be found by matching this limiting solution to another solution that is valid further away from $r=0$.
The method of matched asymptotic expansions is often used to solve this type of problems asymptotically, see e.g. Kevorkian & Cole (1996) "Multiple Scale and Singular Perturbation Methods", Springer. In this method you would need to construct two solutions: "inner" solution, valid as $r\to 0$, (it is the solution derived above) and "outer" solution, valid as $r\to\infty$ (which you took to be a plane wave). Some ingenuity may be needed to ensure that the ranges of validity of these two solutions overlap, so that they can actually be matched.
