I am studying a problem in high energy physics and have reduced it to the following problem regarding the confluent hypergeometric differential equation: $$ \left(\gamma \partial_\gamma^2 + \partial_\gamma - \frac{\ell^2}{4 \gamma} - \frac{\gamma}{4} \right) u(\gamma) = \left(V(\gamma) - E \right) u(\gamma) $$ with $\gamma \in [0,\infty)$ and $V(\gamma)$ any smooth function that interpolates between some value $\Delta_1$ as $\gamma \to 0$ and some other value $\Delta_2$ as $\gamma \to \infty$. We could consider $V(r) = \Delta \tanh(\gamma - \gamma_0)$ as an example (so that $\Delta_1 = 0, \Delta_2 = \Delta$). Alternatively, we could forgo the smoothness assumption and consider a step-function potential such that $V(\gamma) = \Delta_1 \theta(\gamma_0-\gamma) + \Delta_2 \theta(\gamma - \gamma_0)$ for some $\gamma_0 > 0$.
I would like to find the an approximate expression for the eigenvalues $E$ with the boundary conditions specified by requiring a normalizable solution as $\gamma \to 0$ and $\gamma \to \infty$.
Here is the progress I have made so far: for a hard-wall boundary condition i.e., $\Delta_2 = \infty$, the boundary condition becomes $u(\gamma_0) = 0$. In this case, requiring a well-behaved solution near the origin forces one to pick the solution $$u(\gamma) = e^{-\gamma/2} \gamma^{|\ell|/2}\mathcal{M}\left( \frac{1+|\ell|}{2} + \Delta_1 - E,1+|\ell|,\gamma \right)$$ and then the eigenvalues can be found by imposing the boundary condition at $\gamma_0$. Here, I can use the asymptotic expansion of $\mathcal{M}(a,b,z)$ for large $a$ or $z$ and obtain approximate expressions for E in various regimes. (For large $z$, these lead to those of a 2D quantum harmonic oscillator with a parabolic confinement that is placed in a disc geometry). This method gives results that match the expressions derived from the standard WKB type arguments.
However, I am unsure how to proceed in the general case and would appreciate any input. Thus far, what I have attempted is as follows: for $\gamma < \gamma_0$, I pick the solution given by $\mathcal{M}$ as this is regular at the origin. Then, for $\gamma > \gamma_0$ I pick the solution given by $\mathcal{U}$, which is regular at infinity. My guess had been that matching these solutions and their derivatives near $\gamma_0$ would lead to the quantization condition on E, but I have not been able to obtain an approximate expression in this case. We can even set $\Delta_1 = 0$ if this helps, in which case this reduces to the problem of a 2-dimensional quantum oscillator in a magnetic field with a finite step-function potential. I was also unable to locate any literature that addressed this problem and discussed an asymptotic solution.
To match the solutions, what I tried was using the asymptotic expansion for $\mathcal{M}(a,b,z)$ for large $z$ and that of $U(a,b,z)$ for small $z$, and then matching the approximate solutions and their derivatives at $\gamma_0$. I would appreciate any help towards obtaining a solution based on either this approach or even one that uses the semi-classical WKB approximation.
EDIT: As a first step, it would also be helpful to understand how I should find zero-mode solutions $E \approx 0$. From numerical considerations, I expect that these are localized in the vicinity of the boundary $\gamma \approx \gamma_0$. We can also set $\Delta_1 = 0$ in case that simplifies the problem.
EDIT: As suggested in the comments, it would be useful to write the original radial Schrodinger equation of interest. In this case, it is: $$ \hbar^2 \psi''(x) - \left(\hbar^2 \frac{l^2-1/4}{x^2} - \hbar \omega_c \ell + \frac{\omega_c^2 x^2}{4} - 2 (E - V(x))\right) \psi(x) = 0 $$ where again $V(x)$ is a smooth function which interpolates between two values $\Delta_1$ and $\Delta_2$ as $x$ goes from $0$ to $\infty$. Here, $\hbar$ should be identified as the small parameter $\epsilon$.