I am looking at a problem that has a potential
$$ V(r) = \begin{cases} 0 & a<r<b\\ \infty & \text{elsewhere} \end{cases} $$
This is a modification of the infinite spherical well problem, since the well does not extend to $r = 0$. The eigenfunctions are of the form $\Psi(r,\theta,\phi) = R(r)Y_{l,m}(\theta,\phi)$, where the $Y_{l,m}$'s are the usual spherical harmonics. From this, it can be shown that the radial equation becomes
$$ \frac{d^2R}{d\rho^2} + \frac{2}{\rho}\frac{dR}{d\rho} + \bigg[1-\frac{l(l+1)}{\rho^2}\bigg]R = 0, $$ where $\rho\equiv kr$ is a dimensionless parameter, and $k^2 \equiv 2mE/\hbar^2$. This equation is the spherical Bessel equation, and has known solutions that are the spherical Bessel (regular) and spherical Neumann (irregular) functions. The Bessel functions are regular in that they go to zero as $r$ goes to zero, whereas the Neumann functions do not.
This is where I am stuck; can we in this case assume that the solution are just the regular Bessel functions? In the infinite spherical well, where have $V(r) = 0$ for $0<r<a$, so we must have regularity at the origin. In that case the radial functions would be $$ j_l(\rho) = (-\rho)^l \bigg(\frac{1}{\rho}\frac{d}{d\rho}\bigg)^l\frac{\sin\rho}{\rho} $$ So, for example, the ground state function, i.e. $j_0(\rho) = \sin\rho/\rho$. Applying the boundary condition at $r = a$, i.e. $j_0(\rho = ka) = 0$ gives $$ \frac{\sin ka}{ka} = 0 \Rightarrow ka = n\pi $$ Substituting the expression for $k$, we ultimately get $$ E_{n,l=0} = \frac{\pi^2\hbar^2}{2ma^2}n^2 $$
This is all very well known. However, in the problem at hand, we also require that $R(\rho = kb) = 0$ in addition to $R(\rho = ka) = 0$. Thus, I am not sure that we can just use the regular Bessel functions are eigenfunctions. I suppose you can just make the substitution $a \rightarrow b-a$ so that the ground state energies are $$ E_{n,l=0} = \frac{\pi^2\hbar^2}{2m(b-a)^2}n^2, $$ but this is not rigorous at all and seems fishy to me. My question is then: How would you write the radial eigenfunctions in this case, and subsequently how to get the energy levels using both boundary conditions?