Your radial equation is wrongly formulated. It looks as an equation for flat 1D space, lacking a first-derivative term.

Instead of postulating the equation, as you seem to have done, let's derive it. I'll use units such that $\frac{\hbar^2}{2m}=1$, so that the equations look simpler, you should be able to reproduce it with your units.

So, the Schrödinger equation for a particle in a spherically symmetric potential $V(r)$ with $r=|\vec x|$ is:

$$-\nabla^2\Psi(\vec x)+V(r)\Psi(\vec x)=E\Psi(\vec x)$$

As your potential is spherically symmetric, we can make use of this symmetry and switch to spherical coordinates. Our Laplacian would look like:

$$\nabla^2 f=\frac1{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial f}{\partial r}\right)+\frac1{r^2\sin\theta}\frac{\partial}{\partial \theta}\left(\sin\theta\frac{\partial f}{\partial \theta}\right)+\frac1{r^2\sin^2\theta}\frac{\partial^2 f}{\partial\phi^2}$$

With such a Laplacian, our Schrödinger equation is separable, so we can look for solution in the form of $\Psi(r,\theta,\phi)=u(r)v(\theta,\phi)$.

Substituting it into new formula for Laplacian, we have:

$$\nabla^2\Psi(r,\theta,\phi)=v(\theta,\phi)\frac1{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial u(r)}{\partial r}\right)+u(r)\frac1{r^2\sin\theta}\frac{\partial}{\partial \theta}\left(\sin\theta\frac{\partial v(\theta,\phi)}{\partial \theta}\right)+\\+u(r)\frac1{r^2\sin^2\theta}\frac{\partial^2 v(\theta,\phi)}{\partial\phi^2}$$

Or, denoting radial part of Laplacian with $-\hat R$ and angle part with $-\hat L^2/r^2$, we have:

$$\nabla^2\Psi(r,\theta,\phi)=-v(\theta,\phi)\hat R u(r)-u(r)\frac{\hat L^2v(\theta,\phi)}{r^2}$$

Now we can write our Schrödinger equation as:

$$v(\theta,\phi)\hat R u(r)+u(r)\frac{\hat L^2v(\theta,\phi)}{r^2}+V(r)u(r)v(\theta,\phi)=Eu(r)v(\theta,\phi)$$

Multiply both sides by $\frac{r^2}{u(r)v(\theta,\phi)}$ and rearrange the terms:

$$-\frac{r^2\hat Ru(r)}{u(r)}-V(r)+E=\frac{\hat L^2v(\theta,\phi)}{v(\theta,\phi)}$$

Now we have separated radial variables from angle ones, so we introduce a separation constant, which we'll write as $l(l+1)$. It's an eigenvalue of $\hat L^2$ operator (and eigenfunctions are spherical harmonics). Now, multiplying everything by $u(r)/r^2$, we have:

$$\hat Ru(r)+V(r)u(r)+\frac{l(l+1)}{r^2}u(r)=Eu(r),$$

or, finally, writing out expression for $\hat R$,

$$-\frac1{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial u(r)}{\partial r}\right)+V(r)u(r)+\frac{l(l+1)}{r^2}u(r)=Eu(r).$$

This is the equation you should be trying to solve (up to units). Now its solution should be in terms of spherical Bessel functions.

Your radial equation is wrongly formulated. It looks as an equation for flat 1D space, lacking a first-derivative term, arising from spherical coordinates Jacobian.

Instead of postulating the equation, as you seem to have done, let's derive it. I'll use units such that $\frac{\hbar^2}{2m}=1$, so that the equations look simpler, you should be able to reproduce it with your units.

So, the Schrödinger equation for a particle in a spherically symmetric potential $V(r)$ with $r=|\vec x|$ is:

$$-\nabla^2\Psi(\vec x)+V(r)\Psi(\vec x)=E\Psi(\vec x)$$

As your potential is spherically symmetric, we can make use of this symmetry and switch to spherical coordinates. Our Laplacian would look like:

$$\nabla^2 f=\frac1{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial f}{\partial r}\right)+\frac1{r^2\sin\theta}\frac{\partial}{\partial \theta}\left(\sin\theta\frac{\partial f}{\partial \theta}\right)+\frac1{r^2\sin^2\theta}\frac{\partial^2 f}{\partial\phi^2}$$

With such a Laplacian, our Schrödinger equation is separable, so we can look for solution in the form of $\Psi(r,\theta,\phi)=u(r)v(\theta,\phi)$.

Substituting it into new formula for Laplacian, we have:

$$\nabla^2\Psi(r,\theta,\phi)=v(\theta,\phi)\frac1{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial u(r)}{\partial r}\right)+u(r)\frac1{r^2\sin\theta}\frac{\partial}{\partial \theta}\left(\sin\theta\frac{\partial v(\theta,\phi)}{\partial \theta}\right)+\\+u(r)\frac1{r^2\sin^2\theta}\frac{\partial^2 v(\theta,\phi)}{\partial\phi^2}$$

Or, denoting radial part of Laplacian with $-\hat R$ and angle part with $-\hat L^2/r^2$, we have:

$$\nabla^2\Psi(r,\theta,\phi)=-v(\theta,\phi)\hat R u(r)-u(r)\frac{\hat L^2v(\theta,\phi)}{r^2}$$

Now we can write our Schrödinger equation as:

$$v(\theta,\phi)\hat R u(r)+u(r)\frac{\hat L^2v(\theta,\phi)}{r^2}+V(r)u(r)v(\theta,\phi)=Eu(r)v(\theta,\phi)$$

Multiply both sides by $\frac{r^2}{u(r)v(\theta,\phi)}$ and rearrange the terms:

$$-\frac{r^2\hat Ru(r)}{u(r)}-V(r)+E=\frac{\hat L^2v(\theta,\phi)}{v(\theta,\phi)}$$

Now we have separated radial variables from angle ones, so we introduce a separation constant, which we'll write as $l(l+1)$. It's an eigenvalue of $\hat L^2$ operator (and eigenfunctions are spherical harmonics). Now, multiplying everything by $u(r)/r^2$, we have:

$$\hat Ru(r)+V(r)u(r)+\frac{l(l+1)}{r^2}u(r)=Eu(r),$$

or, finally, writing out expression for $\hat R$,

$$-\frac1{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial u(r)}{\partial r}\right)+V(r)u(r)+\frac{l(l+1)}{r^2}u(r)=Eu(r).$$

This is the equation you should be trying to solve (up to units). Now its solution should be in terms of spherical Bessel functions.