The classical harmonic oscillator has energy

$$E_{osc.}=\frac{p^2}{2m}+\frac12m\omega^2x^2\ ,\quad\omega=\sqrt{\frac{k}{m}}$$

When the harmonic oscillator was first studied by physicist under the lens of the quantum realm, there was the need to construct a viable Hamiltonian.

$$H_{osc.}=\frac{\hat p^2}{2m}+\frac12m\omega^2\hat x^2\ ,\quad\left[\hat x,\hat p\right]=i\hbar$$

It's a quantum mechanical system inspire from classical mechanics in the simplest and purest way. Studying it one can see that the oscillator has specific allowed values of amplitude and energy.

Results from such a study can be used almost anywhere. Any system with an equilibrium point has an oscillator-like potential in the neighbouring region, at least in approximation. Meaning that all systems in equilibrium can be approximated up to a degree with a harmonic oscillator.

Just like the infinite well, the potential blows up to infinity, and as such all allowed states of a particle are bound states.

The eigenvalue equation is: $$\hat H\psi_n(x)=E_n\psi_n(x)$$

There are two ways to solve the eigenvalue problem. The first is a brute force solution of the differential equation, and in the process finding out the quantization of energy directly from the form of the equation. The second way is to solve it algebraically by defining the appropriate creation and annihilation operators, which move the system to the next or previous energy state respectively.

The complete equation is:

$$\boxed{-\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2}+\frac12m\omega^2x^2\psi(x)=E\psi(x)} \tag{1}$$

The first obstacle we find when we try to solve it is all those quantities with their respective units. There are two practical ways to get rid of them. The first one is the use of natural units, $\hbar=m=\omega=1$, and then some kind of multiplication between them and the results in order to recreate the proper physical units. Then the problem reduces to the form below:

$$\psi''(x)+\left(2E-x^2\right)\psi(x)=0 \tag{2}$$

The *problem* with this approach when first learning about the quantum harmonic oscillator is that you kind have to already know before what type of result you expect and what unit reconstruction to do, or blindly trust your textbook, ending up in not understanding what you're doing (*that's what happened to me at least*).

The second approach is starting by demanding the use of a dimensionless variable, and then find what combination of constants makes $x$ dimensionless.

We write the variable $x$ as $x=\alpha u$, where $u$ is dimensionless and $\alpha$ has units of length. The kinetic energy part of the equation can be now written as:

$$ -\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2} =-\frac{\hbar^2}{2m}\frac{d^2\psi}{\left(d\left(\alpha u\right)\right)^2} =\frac{\hbar^2}{m\alpha^2}\cdot\left(-\frac12\frac{d^2\psi}{du^2}\right) $$

Similarly, the potential energy part is:

$$ \frac12m\omega^2x^2\psi(x) =\frac12m\omega^2\left(\alpha u\right)^2\psi(u) =m\omega^2\alpha^2\cdot\left(\frac12u^2\psi(u)\right) $$

On the other side of the equation we only have $E\psi(u)$, and hence we have two ways to construct the energy from the constants.

$$ E=\frac{\hbar^2}{m\alpha^2}=m\omega^2\alpha^2 \Leftrightarrow \boxed{\alpha^2=\frac{\hbar}{m\omega}} $$

Putting this back to equation $(1)$ we have: $$- \frac12\frac{\hbar^2}{ma^2}\frac{d^2\psi}{du^2} +\frac12m\omega^2a^2u^2\psi(u) =E\psi(u) $$ $$ -\frac12\hbar\omega\frac{d^2\psi}{du^2} +\frac12\hbar\omega u^2\psi(u) =E\psi(u) $$ $$ -\frac{d^2\psi}{du^2} + u^2\psi(u) = \frac{2E}{\hbar\omega}\psi(u) $$

The left part of the equation is dimensionless. This forces us to interpret the quantity $2E/\hbar\omega$ as also dimensionless. For ease of calculations we consider it as some dimensionless energy $\mathcal E$.

$$ \boxed{\mathcal{E}:=\frac{2E}{\hbar\omega}} $$

Because we included the factor of $2$ in the dimensionless energy, this will result in it missing from equation $(2)$.

$$\boxed{\frac{d^2\psi}{du^2}=\left(u^2-\mathcal{E}\right)\psi(u)}\tag{3}$$

There doesn't seem to be some obvious solution for $(3)$. In reality a solution doesn't exist for all possible energies, since the squared term blows the wavefunction to infinity for large enough $u$, making it impossible to normalize.

We can find the form of the solutions by studying how the equation behaves when $u\to\pm\infty$.

$$\frac{d^2\psi}{du^2}=u^2\psi(u)$$

It seems like there isn't a pure polynomial solution, since the two parts of the equation will have a four degree difference.

One can notice that the exponential $\left(e^{u^2}\right)'=2ue^{u^2}$ can increase the polynomial degree, and then try a solution in the form of $\phi$ below.

$$\phi(u)=u^ke^{\lambda u^2}$$ $$\phi'(u)=2\lambda u^{k+1}e^{\lambda u^2}+\dots = 2\lambda u\cdot\phi(u)+\dots$$ $$\phi''(u)=(2\lambda u)^2\cdot \phi(u)+\dots$$

For $\lambda=\pm1/2$ we have the form we need. So the solution could have the form:

$$\phi(u)\sim Au^{k}e^{-u^2/2}+Bu^ke^{u^2/2}$$

Now the problem of the solution blowing up to infinity is obvious. The second term of the sum cannot exist as a part of the wavefunction due to normalization constraints. We can suppose the form of the solution below:

$$\boxed{\psi(u)=h(u)\cdot e^{-u^2/2}} \tag{4}$$

$\phi(u)$ had a $Au^k$ term where the power $k$ was irrelevant, so we can interpret $h(u)$ as a power series.

$$h(u)=\sum_{n=0}^{\infty}a_nu^n$$

Replacing $(4)$ in $(3)$ will result in a differential equation with respect to $h(u)$, and solving it we can find the wavefuntion $\psi(u)$.

$$\frac{d^2\psi}{du^2}=(u^2-\mathcal{E})\psi(u)$$ $$\left[ \frac{d^2h}{du^2}-2u\frac{dh}{du}+(u^2-1)h(u) \right]e^{-u^2/2} =(u^2-\mathcal{E})h(u)e^{-u^2/2} $$ $$\boxed{ \frac{d^2h}{du^2}-2u\frac{dh}{du}+(\mathcal{E}-1)h=0 } \tag{5} $$

A simple replacement will result in a power series derivative, shifted starting point and a bunch of algebraic manipulation we can avoid. For a random $n\ge0$ some part of the power series will be:

$$h(u)=\ ...+a_nu^n+a_{n+1}u^{n+1}+a_{n+2}u^{n+2}+\dots$$

We can find the parts of the first and second derivative that end up with the same power and then sum them up accordingly.

$$\frac{d^2h}{du^2}=(n+2)(n+1)a_{n+2}u^n$$ $$-2u\frac{dh}{du}=-2na_nu^n$$ $$(\mathcal{E}-1)h=(\mathcal{E}-1)a_nu^n$$

Summing up the three terms above we can solve the DE and find the power series terms.

$$ \sum_{n=0}^\infty\Big[ (n+2)(n+1)a_{n+2}-2na_n+(\mathcal E-1)a_n \Big]\cdot u^n=0 $$ $$(n+2)(n+1)a_{n+2}-(2n+1-\mathcal E)a_n=0$$ $$(n+2)(n+1)a_{n+2}=(2n+1-\mathcal E)a_n$$

$$\boxed{a_{n+2}=\frac{2n+1-\mathcal E}{(n+2)(n+1)}a_n} \tag{6}$$

With this recurrence relation we can find the whole function $h(u)$ if we only knew the first two terms, $a_0$ and $a_1$. The first one will give us all the even part of $h(u)$, and the second the odd part. In order to find those terms, zeroing the function and its first derivative is enough.

$$h(0)=a_0+a_1\cdot0+a_2\cdot0^2+\ ...\ =a_0$$ $$h'(0)=a_1+2a_2\cdot0+3a_3\cdot0^2+\ ...\ =a_1$$

There still exists a problem. The series doesn't always converge. A simple ratio test can show that:

$$\left.\frac{a_{n+2}}{a_n}\right|_{n\to+\infty}=\frac{2n}{n^2}=\frac2n$$

The terms do get smaller, but not fast enough for convergence. It might not be obvious initially, but it can be shown that this is similar scaling to that of $e^{u^2}$, which for obvious reasons diverges.

$$\begin{align*} e^{u^2} &= \sum_{n=0}^\infty\frac{1}{n!}\left(u^2\right)^n \\ &= \sum_{n=0}^\infty\frac{1}{n!}u^{2n} \\ {\small(j=2n)}&=\sum_{j=0,2,4,\small\dots}^\infty\frac{1}{(j/2)!}u^j \\ &= \sum_{j=0,2,4,\small\dots}^\infty a_ju^j,\quad a_j=\frac{1}{(j/2)!} \end{align*}$$ We have to keep in mind since $j$ is even, the factorial of $j/2$ is the factorial of a simple integer, so nothing fishy is happening here. Now, studying the limit at infinity we that we have the same scaling behaviour as above.

$$ \frac{a_{j+2}}{a_j} =\frac{(j/2)!}{\big((j+2)/2\big)!} =\frac{(j/2)!}{(j/2)!\cdot\big((j+2)/2\big)} =\frac{1}{\big((j+2)/2\big)} =\frac{2}{j+2} $$ $$\left.\frac{a_{j+2}}{a_j}\right|_{j\to+\infty}=\frac{2}{j}$$

We just proved that if the power series has infinite terms, the wavefunction at infinities would behave like below:

$$h(u)\propto e^{u^2},$$ $$\psi(u)\propto e^{u^2}\cdot e^{-u^2/2}=e^{u^2/2},\quad u\to\pm\infty \tag{7}$$

In such case the wavefunction can't be normalized, and thus cannot exist physically. The only way $h(u)$ has an acceptable form is if the power series terminates at some point, meaning $2n+1-\mathcal E=0$. Thus, the energy is quantized.

$$\boxed{\mathcal E_n=2n+1} \tag{8}$$

Replacing our definition for the dimensionless energy $(\mathcal E:=2E/\hbar\omega)$ we find the actuall energy levels in the real world.

$$\boxed{E_n=\left(n+\frac12\right)\hbar\omega} \tag{9}$$