I had a question on the "correct" numerical discretization of the Laplacian for differential eigenvalue problems on the real line. As a simple model for the kind of problem I have consider a non-self-adjoint Schrödinger operator eigenvalue problem like $$-u_{xx} + u_x + q(x) u = \lambda u$$ on the whole line: $u \in L_2({\mathbb R}).$ (Here $q(x)$ is as smooth and rapidly decaying as you like -- say Schwartz class). My actual problem is somewhat more complicated but I think that the above captures all of the important features. I would like to compute a numerical approximation to the spectrum of this operator. It is clear from a Weyl sequence argument that the essential spectrum for the above problem is $\sigma_{\text{ess}}=\{\lambda ~\vert~ \lambda = i k + k^2 ~~~~k\in {\mathbb R} \}.$
My naive approach to computing the spectrum numerically is to just pose the problem on some large but finite interval $[-L,L]$ and to discretize the problem on that interval. The thing that I quickly realized is that form of the discretization can dramatically change the result. For instance the naive discretization would be to replace the second derivative with the tridiagonal matrix $(\Delta x)^{-2} {\bf M}_{(2)}$, where the matrix ${\bf M}_{(2)}$ is tridiagonal matrix with $2$ along the diagonal and $-1$ the sub and sup diagonals, and to replace the first derivative with the matrix $(2\Delta x)^{-1}{\bf M}_{(1)}$, where ${\bf M}_{(1)}$ has zero along the main diagonal, $1$ along the sup diagonal and $-1$ along the sub diagonal. Unfortunately this discretization (corresponding to Dirichlet boundary conditions at the endpoints of the interval) gets the essential spectrum completely wrong: it is easy to see that the spectrum of ${\bf M}_{(2)}+ {\bf M}_{(1)}$ is purely real as long as $\Delta x$ is less than $2$. Neumann boundary conditions suffer from the same problem, but the periodic versions of the first and second derivatives get the essential spectrum correct. I had naively assumed that on a large box the different boundary conditions would give similar spectra, but this is very much not true.
I think that I more or less understand this particular example but my question is whether there is some kind of theory in numerical analysis for the "correct" discretization of differential operators like the Laplacian on the whole line. The whole question seems to be closely related to the question of different self-adjoint extensions of an unbounded operator and deficiency indices a la Reed and Simon.