We are given Hermitian operator of the form
$H(x) =(-\hbar^2/2m) \partial^2/\partial x^2 + V(x)$
(where $\hbar$ and $m$ are real constants) which has orthogonal eigenfunctions corresponding to a set of distinct real eigenvalues, whose size is finite or countably infinite.
We also are given
$-i\hbar\partial\Psi(x,t)/\partial t = H(x)\Psi(x,t)$
(Schrodinger’s Equation, of quantum mechanics, where $i = \sqrt{-1}$)
We wish to solve for $\Psi$. Also, $\Psi$ is restricted to be a “test function”, in the distribution theory sense.
I have been told that, if $\Psi$ is a solution to the above equation, it must lie in the function space spanned by a linear combination of the eigenfunctions of $H$. Is that mathematically provable to be always true, sometimes true, or maybe just a postulate of QM so that, although it is not mathematically always true, we assume it is always true for physical systems?
I have not looked into the case of $\mathbb{R}$, but I have looked into the case on $[0,\infty)$ (or $(-\infty,0]$). In the case on $[0,\infty)$, for example, the Sturm-Liouville operator $H$ with no singular points in $(0,\infty)$, where you have a possibe boundary condition at $0$ of the form $$ \cos\alpha \psi(0)+\sin\alpha \psi'(0)=0, $$ you can solve for the unique classical eigenfunctions $\psi_{\lambda}$ satisfying \begin{align} \cos\alpha \psi_{\lambda}(0)+\sin\alpha \psi_{\lambda}'(0)&=0, \\ -\sin\alpha \psi_{\lambda}(0)+\cos\alpha \psi_{\lambda}(0)&=1. \end{align} (These are classical ODE solutions that may lie in $L^2[0,\infty)$ for some $\lambda$, but generally do not unless $\lambda$ is in the discrete spectrum.) When you do this, you find that there is a unique positive spectral density measure $\rho$ such that the following outer integral converges in $L^2[0,\infty)$: $$ f = \int_{-\infty}^{\infty} \left(\int_{0}^{\infty}f(x')\psi_{\lambda}(x')dx'\right) \psi_{\lambda}(x) d\rho(\lambda). $$ There is a corresponding Plancherel identity, $$ \|f\|^2_{L^2[0,\infty)} = \int_{-\infty}^{\infty}\left|\int_{0}^{\infty}f(x')\varphi_{\lambda}(x')dx'\right|^2 d\rho(\lambda). $$ In this setting it is true that the every function in $L^2[0,\infty)$ can be built up from the eigenfunctions of the Hamiltonian. This case fits rather nicely with Dirac's formulation; it fact, Dirac probably borrowed much of his formalism from the results obtained for Sturm-Liouville Theory in the early part of the 20th century by Weyl and others.
I have not thoroughly looked into how to link the two results on each half-line, but it is my understanding this can be done (possibly with a $2\times 2$ matrix?). So I would say there is a sense in which this statement of yours is true: "I have been told that, if Ψ is a solution to the above equation, it must lie in the function space spanned by a linear combination of the eigenfunctions of H."