Let $H$ be a Schrödinger operator on $L^2(\mathbb R)$, that is a self-adjoint operator, such that $$ (Hf)(x) = -f''(x) + V(x) \, f(x) \qquad \text{for } f \in \mathrm D(H) \: . $$ These operators play a huge role in Quantum Mechanics – especially their spectrum, which represents the possible energies a system can have. To find the pure point spectrum (ie. eigenvalues), one simply needs to solve the equation $$ Hf = \epsilon f, \qquad \text{where} \quad f \in \mathrm D(H), \quad \epsilon \in \mathbb R \: , \tag{TISE} \label{TISE} $$ in physics called the Time-independent Schrödinger Equation. However, to find the continuous spectrum, one has to use more advanced techniques from Functional Analysis. That being said, physicists quite frequently “find” the spectrum by solving \eqref{TISE}, but taking $f$ which are not square integrable. For example the TISE of a free Hamiltonian ($V=0$) is: $$ -f'' - \epsilon f = 0 \qquad \implies \qquad f(x) = A \, e^{-\mathrm i \, \sqrt{\epsilon} \, x} + B \, e^{+\mathrm i \, \sqrt{\epsilon} \, x} \quad \text{for} \quad \epsilon \in [0, \infty) \: . $$ Peculiarly, the spectrum of the free Hamiltonian really is $[0, \infty)$. A similar things happens in the derivation of Bloch's theorem, when $f$ is taken to be a periodic function (up to a phase) – such a wavefunction is definitely not square-integrable, but it leads to a correct prediction of the spectrum of $H$.
This doesn't seem to be just a coincidence. Is there a (mathematically rigorous) theorem, which would guarantee that the solutions of \eqref{TISE} determine the spectrum of $H$, even though they are not in the domain of $H$? When I ask the physicists, they usually respond “[something something] Rigged Hilbert space.” but I couldn't find any general result which would justify this approach.
Any valid Hamiltonian is self-adjoint and has a basis of eigenfunctions, which may include discrete and continuous spectral components. You can expand functions in $L^2(\mathbb{R})$ with respect to the eigenfunctions associated with some other Hamiltonian because the basis will be complete. It comes down to choosing a convenient set of eigenfunctions for the expansion. The reference coordinates are spectral coordinates for some canonical problem; the eigenfunctions for the canonical problem are sufficient to expand anything else, but the problem should be reduced by choosing a canonical type of problem and its associated eigenfunctions.