I want to solve the following eigenvalue problem
$$ \left(-\frac{\partial^2 }{\partial x^2 } + x^2 \right) \psi = E \psi , $$
with the boundary condition $ \psi'(0) =\beta \psi(0 ) $, and $\psi(x\rightarrow \infty) = 0 $. Here $\beta $ is an arbitrary real constant.
It is well known that the spectral method is much more accurate than the common finite-difference method, so I want to adopt the former in this problem. As far as I understand, the essence of the spectral method is to choose a basis (not necessarily orthogonal) for the Hilbert space and then construct the hamiltonian matrix $H$ and the overlap matrix $ S$ in this basis and then solve the generalized eigenvalue problem $Hv = E S v $.
The question is, how to choose or construct a basis for the problem?
Theoretically, any basis set of the underlying Hilbert space should work. As noted in the comments, it makes sense to take your domain into considerations, since different basis sets are defined on different domains.
Convergence patterns, i.e., how many functions you need from the underlying basis to have a good approximation of your eigenfunctions, highly depend on your choice of the basis. To achieve fast convergence you should use a basis set of functions that solves a similar problem. Note that your operator, i.e., $H = \frac{\partial ^2}{\partial x^2}+x^2$, is very similar to the operator of the quantum harmonic oscillator (QHO) problem. Hence, it makes sense to use as a basis Hermite functions, which are the eigenfunctions of the QHO.