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.
