2D isotropic quantum harmonic oscillator: polar coordinates This may come a bit elemental, what I was working on a direct way to find the eigenfunctions and eigenvalues of the isotropic two-dimensional quantum harmonic oscillator but using polar coordinates:
$$
H=-\frac{\hbar}{2M}\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\right)+\frac{M\omega^2}{2}\left(x^2+y^2\right).
$$
I can easily solve the 2-dimensional case in cartesian coordinates as we can separate the hamiltonian in independent oscillators for each coordinate. For the polar case in two dimensions, we can rewrite,
$$H=-\frac{\hbar}{2M}\left(\frac{\partial^2}{\partial r^2}+\frac{1}{r}\frac{\partial}{\partial r}+\frac{1}{r^2}\frac{\partial^2}{\partial \phi^2}\right)+\frac{\omega^2}{2}r^2.$$
With $r^2=x^2+y^2$ and $\phi=\arctan(y/x)$. 
Using separation of variables $\psi(r,\phi)=R(r)\Phi(\phi)$ and plugging into the Schrödinger equation, we can easily solve for the angular part $\Phi=e^{im\phi}$, where $m\in \mathbb{Z}$.
Plugging back into the Schrodinger equation, for the radial part, we get:
$$r^2R''+rR'+(r^2E-m^2-M\omega^2r^4)R=0.$$
While I have an idea for the solution by making an analogy with the 3D case (where we get Laguerre polynomials), I'm not sure how to correctly proceed from here. I appreciate any input or even useful references* (all of the references I've found deal with the 3D case, which I have no problem solving).
*I read online this problem is treated in the book of "Wave Mechanics" from Pauli, but unfortunately it isn't available on neither of my campus libraries nor online (it's only available for purchase and I lack the funds to buy it). 
 A: 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).
A: I think you may have a factor of 2 error for the $E$ and the power of $M$, here is my derivation, I use $\theta$ instead of $\phi$.(we may double check)
In polar coordinates, the Del operator $\nabla^2$ is defined as:
$$
\begin{align}
\nabla^2
&=\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial f}{\partial r}\right)+\frac{1}{r^{2}} \frac{\partial^{2} f}{\partial \theta^{2}}\\
&=\frac{\partial^2 f}{\partial r^2}+\frac{1}{r} \frac{\partial f}{\partial r}+\frac{1}{r^{2}} \frac{\partial^{2} f}{\partial \theta^{2}}
\end{align}
$$
Then the Shrodinger Equation for this system can be written as:
$$
\left(-\frac{\hbar^2\nabla^2}{2M}+\frac{M\omega^2r^2}{2}\right)\Psi(r,\theta)=E\Psi(r,\theta)\\
\left(-\frac{\partial^2}{2M\partial r^2}-\frac{1}{2Mr} \frac{\partial}{\partial r}-\frac{1}{2Mr^2} \frac{\partial^{2} }{\partial \theta^{2}}+\frac{M\omega^2r^2}{2}\right)\Psi(r,\theta)=E\Psi(r,\theta)\\
\left(-\frac{\partial^2}{\partial r^2}-\frac{1}{r} \frac{\partial}{\partial r}-\frac{1}{r^2} \frac{\partial^{2} }{\partial \theta^{2}}+M^2\omega^2r^2\right)\Psi(r,\theta)=2ME\Psi(r,\theta)
$$
By assuming that the solution is separatable $\Psi(r,\theta)=R(r)\psi(\theta)$, we can solve the angular part fairly easily:
$$
\psi(\theta)=e^{im\theta}\qquad m=0,1,2\dots
$$
The the radia part can be rearraged when substitude the angular solution into the Scrodinger equation:
$$
r^2R''+rR'+ \left(2r^2ME-m^2-M^2\omega^2r^4\right)R=0
$$
The equation can be simplfied by setting $M=\omega=1$
$$
r^2R''+rR'+ \left(2r^2E-m^2-r^4\right)R=0
$$
