The hypothetical balls are part of a single quantum system, i.e., there can be (and indeed are) quantum mechanical correlations between them.
If the system is in a state representing a single particle, then it is known that only one ball is excited, but it is uncertain which ball it is.
For each ball, there is a probability amplitude that it is the one that is excited. If you write a function for the probability amplitude that the ball at a particular position is excited, that gives you the quantum wavefunction of the particle.
Indeed, as suggested by phase-space quantization, most of these equations are reducible to generalized Laguerre's, the cousins of Hermite. As universally customary, I absorb $\hbar$, M and ω into r,E. Note your E is twice the energy.
Since $r\geq 0$ you don't lose negative values, and you may may redefine $r^2\equiv x$, so that
$$
r\partial_r = 2x \partial_x \qquad \Longrightarrow r\partial_r (r\partial_r)= r^2\partial_r^2+ r\partial_r=4(x^2\partial_x^2+x\partial_x),
$$
hence your radial equation reduces to
$$
\left ( \partial_x^2+ \frac{1}{x}\partial_x +\frac{E-x}{4x} -\frac{m^2}{4x^2} \right ) R(m,E)=0 ~.
$$
Now, further define
$$
R(m,E)\equiv x^{|m|/2} e^{-x/2} ~ \rho(m,E),
$$
to get
$$
\partial_x R(m,E)= x^{|m|/2} e^{-x/2} \left (-1/2 +\frac{|m|}{2x} + \partial_x \right )~ \rho(m,E) \\
\partial_x^2 R(m,E)= x^{|m|/2} e^{-x/2} \left (-1/2 +\frac{|m|}{2x} + \partial_x \right )^2~ \rho(m,E),
$$
whence the generalized Laguerre equation for non-negative m=|m|,
$$
x \partial_x^2\rho(m,E) +\left({m+1} -x\right )\partial_x \rho(m,E)+\frac{1}{2}(E/2-m-1) \rho(m,E)=0~.
$$
This equation has well-behaved solutions for non-negative integer $$k=(E/2-m-1)/2\geq 0 ~,$$ to wit, generalized Laguerre (Sonine) polynomials $L^{(m)}_k (x)=x^{-m}(\partial_x -1)^k x^{k+m}/k!$.
Plugging into the factorized solution and the above substitutions nets your eigen-wavefunctions. The ground state is $k=0=m$, ($E=2$ in your conventions), so a radially symmetric Gaussian, $e^{-r^2/2}$.
Again, in your idiosyncratic convention, the degeneracy is E/2.
So, degeneracy 2 for $E=4$ : $m=1$, $k=0$; you may check this is just $r e^{-r^2/2 +i\phi} $. You may choose the $\cos \phi$ and $\sin \phi$ solutions, if you wish, constituting a doublet of the underlying degeneracy group SU(2).
Best Answer
Combining Eq. (2.58) and (2.62) gives
\begin{equation*} \hbar \omega \left( a_{\pm } a_{\mp } \pm \frac{1}{2}\right) \psi _{n} =\hbar \omega \left( n+\frac{1}{2}\right) \psi _{n} \end{equation*}
\begin{equation*} \left( a_{\pm } a_{\mp } \pm \frac{1}{2}\right) \psi _{n} =\left( n+\frac{1}{2}\right) \psi _{n} . \end{equation*}
This is two equations in one, one of which is
\begin{equation*} \left( a_{+} a_{-} +\frac{1}{2}\right) \psi _{n} =\left( n+\frac{1}{2}\right) \psi _{n} \end{equation*}
\begin{equation*} \rightarrow a_{+} a_{-} \psi _{n} =n\psi _{n} . \end{equation*}
The other one is
\begin{equation*} \left( a_{-} a_{+} -\frac{1}{2}\right) \psi _{n} =\left( n+\frac{1}{2}\right) \psi _{n} \end{equation*}
\begin{equation*} \rightarrow a_{-} a_{+} \psi _{n} =( n+1) \psi _{n} . \end{equation*}