The latter operator in quantum mechanics are the operators 

$$a^\dagger \ = \  \frac{1}{\sqrt{2}} \left(-\frac{d}{dq} + q\right)$$
and 
$$a \ \ = \  \frac{1}{\sqrt{2}} \left(\ \ \ \!\frac{d}{dq} + q\right).$$

They are differential operators on $\mathbb R.$ If one writes them in the Hermite basis, then 

$$a^\dagger = \begin{pmatrix}
0 & 0 & 0 & 0 & \dots & 0 & \dots \\
\sqrt{1} & 0 & 0 & 0 & \dots & 0 & \dots \\
0 & \sqrt{2} & 0 & 0 & \dots & 0 & \dots \\
0 & 0 & \sqrt{3} & 0 & \dots & 0 & \dots \\
\vdots & \vdots & \vdots & \ddots & \ddots & \dots & \dots \\
0 & 0 & 0 & \dots & \sqrt{n} & 0 & \dots & \\
\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \ddots \end{pmatrix}$$ and 

$$a =\begin{pmatrix}
0 & \sqrt{1} & 0 & 0 & \dots & 0 & \dots \\
0 & 0 & \sqrt{2} & 0 & \dots & 0 & \dots \\
0 & 0 & 0 & \sqrt{3} & \dots & 0 & \dots \\
0 & 0 & 0 & 0 & \ddots & \vdots & \dots \\
\vdots & \vdots & \vdots & \vdots & \ddots & \sqrt{n} & \dots \\
0 & 0 & 0 & 0 & \dots & 0 & \ddots \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \end{pmatrix}.$$

Now assume I was interested in numerically computing the spectrum of 

$$H = \begin{pmatrix} 0& a\\a^* & 0\end{pmatrix}.$$

I absolutely know that this can be computed by hand, but I wonder about how to do this numerically.

A naive idea would be to truncate the above matrices at a large size $N$, but this leads to the wrong spectrum as both matrices then have a non-zero nullspace once they are truncated (it is clear since 0 is then an eigenvalue of geometric multiplicity $1$ for both matrices). Hence, the truncated numerics would predict that the Hamiltonian $H$ has an eigenvalue $0$ of multiplicity 2 rather than 1, which is correct. 

Does anybody know how to numerically overcome this pseudospectral effect?