This is literally the mandatory big matrix problem. We must first find the general expressions for the matrix elements of both the position and the momentum operators for the
harmonic oscillator. Now, from previous problems, we know that we can expression the position and momentum operators in terms of the raising and lowering operators:
$$x \ = \ \sqrt{\frac{\hbar}{2m\omega}} \ (a_{+} \ + \ a_{-}) \ \ \ \text{and} \ \ \ p \ = \ i\sqrt{\frac{\hbar m \omega}{2}} \ (a_{+} \ - \ a_{-})$$
To calculate the matrix elements, we calculate the quantities $\langle n | x | n' \rangle$ and $\langle n | p | n' \rangle$. Let us start with the position operator. We have:
$$\langle n | x | n' \rangle \ = \ \sqrt{\frac{\hbar}{2m\omega}} \ \langle n | (a_{+} \ + \ a_{-}) | n' \rangle \ = \ \sqrt{\frac{\hbar}{2m\omega}} \ (\langle n | a_{+} | n' \rangle \ + \ \langle n | a_{-} | n' \rangle)$$
Also recall from previous questions that:
$$\displaystyle\int f^{*} a_{\pm} g \ \text{dx} \ = \ \displaystyle\int (a_{\mp}f)^{*} g \ \text{dx}$$
This means that we have:
$$\sqrt{\frac{\hbar}{2m\omega}} \ (\langle n | a_{+} | n' \rangle \ + \ \langle n | a_{-} | n' \rangle) \ = \ \sqrt{\frac{\hbar}{2m\omega}} \ (a_{-}\langle n | n' \rangle \ + \ \langle n | a_{-} | n' \rangle) \ = \
\sqrt{\frac{\hbar}{2m\omega}} \ (\sqrt{n}\langle n \ - \ 1 | n' \rangle \ + \ \sqrt{n'}\langle n | n' \ - \ 1 \rangle)$$
$$\Rightarrow \ \sqrt{\frac{\hbar}{2m\omega}} \ (\sqrt{n}\langle n \ - \ 1 | n' \rangle \ + \ \sqrt{n'}\langle n | n' \ - \ 1 \rangle) \ = \
\sqrt{\frac{\hbar}{2m\omega}} \ (\sqrt{n}\delta_{n \ - \ 1, \ n'} \ + \ \sqrt{n'}\delta_{n, \ n' \ - \ 1})$$
We follow a similar procedure to get the matrix elements of the momentum operator:
$$\langle n | p | n' \rangle \ = \ i\sqrt{\frac{\hbar m \omega}{2}} \ \langle n | (a_{+} \ - \ a_{-}) | n' \rangle \ = \ i\sqrt{\frac{\hbar m \omega}{2}} \ (\langle n | a_{+} | n' \rangle \ - \ \langle n | a_{-} | n' \rangle)
\ = \ i\sqrt{\frac{\hbar m \omega}{2}} \ (\sqrt{n}\delta_{n \ - \ 1, \ n'} \ - \ \sqrt{n'}\delta_{n, \ n' \ - \ 1})$$
Now, I'm not sure if this next part constitutes a "completely formal proof", but I think it should be Ok for the purposes of this question. From the matrix elements, we can see that
the non-zero values of both matrices will be directly "off" the main diagonal. We must calculate the Hamiltonian, which is given by:
$$H \ = \ \frac{1}{2m} P^2 \ + \ \frac{m\omega^2}{2} X^2$$
For the matrix $X$, from the general formula for matrix elements, we have:
$$X \ = \ \sqrt{\frac{\hbar}{2m\omega}} \ \begin{pmatrix}
0 & \sqrt{1} & 0 & 0 & \dots \\
\sqrt{1} & 0 & \sqrt{2} & 0 & \dots \\
0 & \sqrt{2} & 0 & \sqrt{3} & \dots \\
0 & 0 & \sqrt{3} & 0 \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{pmatrix}$$
$$X^2 \ = \ \frac{\hbar}{2m\omega} \ \begin{pmatrix}
0 & \sqrt{1} & 0 & 0 & \dots \\
\sqrt{1} & 0 & \sqrt{2} & 0 & \dots \\
0 & \sqrt{2} & 0 & \sqrt{3} & \dots \\
0 & 0 & \sqrt{3} & 0 \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{pmatrix} \begin{pmatrix}
0 & \sqrt{1} & 0 & 0 & \dots \\
\sqrt{1} & 0 & \sqrt{2} & 0 & \dots \\
0 & \sqrt{2} & 0 & \sqrt{3} & \dots \\
0 & 0 & \sqrt{3} & 0 \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{pmatrix} \ = \
\frac{\hbar}{2m\omega}\begin{pmatrix}
1 & 0 & \sqrt{2} & 0 & \dots \\
0 & 3 & 0 & \sqrt{6} & \dots \\
\sqrt{2} & 0 & 5 & 0 & \dots \\
0 & \sqrt{6} & 0 & 7 \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{pmatrix}$$
Notice that for each sum on the diagonal, we are summing two consecutive numbers. Since this matrix is infinite in size, this sum of the two consecutive numbers
will continue on indefinitely. For the matrix $P$, also using the general matrix elements, we have:
$$P \ = \ i\sqrt{\frac{\hbar m \omega}{2}} \begin{pmatrix}
0 & -\sqrt{1} & 0 & 0 & \dots \\
\sqrt{1} & 0 & -\sqrt{2} & 0 & \dots \\
0 & \sqrt{2} & 0 & -\sqrt{3} & \dots \\
0 & 0 & \sqrt{3} & 0 \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{pmatrix}$$
$$P^2 \ = \ -\frac{\hbar m \omega}{2} \begin{pmatrix}
0 & -\sqrt{1} & 0 & 0 & \dots \\
\sqrt{1} & 0 & -\sqrt{2} & 0 & \dots \\
0 & \sqrt{2} & 0 & -\sqrt{3} & \dots \\
0 & 0 & \sqrt{3} & 0 \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{pmatrix}\begin{pmatrix}
0 & -\sqrt{1} & 0 & 0 & \dots \\
\sqrt{1} & 0 & -\sqrt{2} & 0 & \dots \\
0 & \sqrt{2} & 0 & -\sqrt{3} & \dots \\
0 & 0 & \sqrt{3} & 0 \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{pmatrix} \ = \ \frac{\hbar m \omega}{2} \begin{pmatrix}
1 & 0 & -\sqrt{2} & 0 & \dots \\
0 & 3 & 0 & -\sqrt{6} & \dots \\
-\sqrt{2} & 0 & 5 & 0 & \dots \\
0 & -\sqrt{6} & 0 & 7 \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{pmatrix}$$
Now we have:
$$H \ = \ \frac{1}{2m} \Big( \frac{\hbar m \omega}{2} \Big) \ \begin{pmatrix}
1 & 0 & -\sqrt{2} & 0 & \dots \\
0 & 3 & 0 & -\sqrt{6} & \dots \\
-\sqrt{2} & 0 & 5 & 0 & \dots \\
0 & -\sqrt{6} & 0 & 7 \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{pmatrix} \ + \ \frac{m\omega^2}{2} \Big( \frac{\hbar}{2m\omega} \Big) \ \begin{pmatrix}
1 & 0 & \sqrt{2} & 0 & \dots \\
0 & 3 & 0 & \sqrt{6} & \dots \\
\sqrt{2} & 0 & 5 & 0 & \dots \\
0 & \sqrt{6} & 0 & 7 \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{pmatrix} \ = \ \frac{\hbar \omega}{4} \ \begin{pmatrix}
2 & 0 & 0 & 0 & \dots \\
0 & 6 & 0 & 0 & \dots \\
0 & 0 & 10 & 0 & \dots \\
0 & 0 & 0 & 14 \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{pmatrix}$$
$$\Rightarrow \ H \ = \ \begin{pmatrix}
\frac{\hbar \omega}{2} & 0 & 0 & 0 & \dots \\
0 & \frac{3\hbar \omega}{2} & 0 & 0 & \dots \\
0 & 0 & \frac{5\hbar \omega}{2} & 0 & \dots \\
0 & 0 & 0 & \frac{7\hbar \omega}{2} \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{pmatrix}$$
And so the general element of the diagonal of the Hamiltonian matrix will be:
$$H_{n} \ = \ \Big(n \ + \ \frac{1}{2} \Big) \hbar \omega$$
Which is exactly what is expected for the Hamiltonian for the quantum harmonic oscillator (since the diagonal is composed of the allowed energies), and each of the
different time-independent solutions of the harmonic oscillator are orthonormal state vectors in $\mathbb{R}^{\omega}$.