Solution: The Schrödinger equation for this system is: $$-\frac{\hbar^2}{2m}\left ( \frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2} \right )\psi(x,y)+\frac{1}{2}m\omega^2(x^2+y^2)\psi(x,y)=E\psi(x,y)$$
We can separate the variables and assume that the wave function can be written as a product of two functions, $\psi(x,y)=X(x)Y(y)$.Substituting this into the Schrödinger equation and dividing by $XY$, we get:
$$-\frac{1}{2m}\frac{X''}{X}=\frac{1}{2m}\frac{Y''}{Y}+\frac{1}{2}m\omega^2(x^2+y^2)=E$$
The LHS of this equation depends only on x, while the RHS depends only on y. Therefore, both sides must be equal to a constant, which we can denote by $E_n$,
$$-\frac{1}{2m}\frac{X''}{X}=E_n , \frac{1}{2m}\frac{Y''}{Y}=E-E_n.$$
These are two independent one-dimensional harmonic oscillator problems, and their solutions are well-known. The energy eigenvalues for the motion in x-direction are:
$$E_n=\hbar\omega\left(n+\frac{1}{2}\right), n=0,1,2,...$$
Similarly, the energy eigenvalues for the motion in y-direction are given by the same formula. Therefore, the total energy eigenvalues for the two-dimensional system are:
$$E_{n_x,n_y}=\hbar\omega\left(n_x+n_y+1\right), n_x,n_y=0,1,2,...$$
The corresponding eigenfunctions are products of the one-dimensional harmonic oscillator wave functions:
$$\psi_{n_x,n_y}(x,y)=\phi_{n_x}(x)\phi_{n_y}(y),$$
where
$$\phi_n(x)= \frac{1}{\sqrt{2^n n!}}\left(\frac{m\omega}{\pi\hbar}\right)^{1/4}H_n \left(\sqrt{\frac{m\omega}{\hbar}}x \right) e^{-m\omega x^2/2\hbar},$$
and $H_n$ are the Hermite polynomials.
