As $r$ approaches zero, the coefficient of $\psi$ becomes dominated by the contribution of $-ar^{-b}$. -a/r^b$. This means that in the vicinity of zero your solution is dominated by the 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). 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 (complex) constant at $r\to0$ depends only on $C_1$ and C_1$, unless $C_2$. C_1$ vanishes. It seems pretty clear that there are no singularities here; typically you should expect that $\psi\sim r$ for $ \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$0. The specific $$ Specific values of $C_1$ and $C_2$ are less trivial to obtain, but they can be found by matching this limiting solution to the 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$. r\to\infty$ (which you took to be a plane wave). Some ingenuity may be needed to ensure that the ranges of vaidity validity of these two solutions overlap, so that they can actually be matchedtogether.

As $r$ approaches zero, the coefficient of $\psi$ becomes dominated by the contribution of $-ar^{-b}$. This means that in the vicinity of zero your solution is dominated by the 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 (complex) constant at $r\to0$ depends on $C_1$ and $C_2$. It seems pretty clear that there are no singularities here; $\psi\sim r$ for $r\to 0$. The specific values of $C_1$ and $C_2$ can be found by matching this limiting solution to the 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$, and "outer" solution, valid as $r\to\infty$. Some ingenuity may be needed to ensure that ranges of vaidity of these two solutions overlap, so that they can be matched together.