Consider the following symmetric (Hermitian) operator on $L^2 (\mathbb{R})$:
$$\hat O= \hat x^3 \hat p +\hat p \hat x^3 \rightarrow \hat O (x)= -i x^3 \partial_x (\cdot) - i\partial_x (x^3 \cdot)$$
It can be proven that the function $$ \psi(x)=\frac{1}{x^{\frac32}}e^{\frac{-1}{4x^2}}$$ satisfies the following eigenvalue equation: $$\hat O \psi(x)=-i \psi(x)$$
Let's see this in a couple of ways: Let us compute the "eigenfunction" of this operator with eigenvalue $-i$:
$$-i x^3 \partial_x \psi -i\partial_x(x^3\psi)=-i\psi$$
$$ x^3 \psi' +3x^2\psi +x^3\psi'=\psi$$
$$ 2x^3 \psi' =(1-3x^2)\psi$$
$$ \frac{\psi'}{\psi} =\frac{1-3x^2}{2x^3}$$
$$\log \psi=-\frac{1}{4x^2}-\frac32\log(x)$$
$$ \psi(x)=\frac{1}{x^{\frac32}}e^{\frac{-1}{4x^2}}$$
And if you don't believe the solution of the differential equation for some reason, by inspection, one can quickly check that $\psi(x)$ satisfies $$\hat O \psi(x)=-i \psi(x)$$
Now, the fact of the matter: Since the operator is Hermitian and this "eigenvalue" is complex, the operator is not self-adjoint, and then $\lambda=-i$ should be in the residual spectrum of $\hat O$.
The problem comes when one checks that the complex-valued function $$ \psi(x)=\frac{1}{x^{\frac32}}e^{\frac{-1}{4x^2}}$$ is smooth and has indeed a finite $L^2$ norm. Using the usual determination (principal argument) for the square root, one quickly gets that $$ |\psi(x)|^2=\frac{1}{|x|^3}e^{\frac{-1}{2x^2}}$$ and the integral on the real line of this yields a finite value (equal to 2).
Now here's the question: The "eigenfunction" corresponding to an element of the residual spectrum should not be in the Hilbert space so I am a bit confused as to why this function seems to be $L^2$.
Can anyone shed some light in this? Any help will be much appreciated :).
Yes, your operator is "symmetric" or "Hermitian" on test functions or on Schwartz functions, because integration by parts succeeds with no boundary term(s). However, for your function, this symmetry just barely fails, because the boundary terms are not $0$, because your function does not decay very well.
Yes, this is a tangible manifestation of the issue of needing to carefully specify domains of unbounded operators...
EDIT: in further detail: for eigenvalue $\lambda$, solutions to the eigenvector equation are linear combinations of $u=e^{i\lambda/4x^2}/|x|^{3/2}$ in regimes $x<0$ and $x>0$. These are in $L^2$ only for $\Im(\lambda)>0$, in which case the solution space is $2$-dimensional (since we can splice different multiples together in $x<0$ and $x>0$). For $\Im(\lambda)\le 0$, there are no $L^2$ solutions.
Since integration by parts works fine for $\int_{\mathbb R} u\cdot \varphi$ for $\varphi$ a Schwartz function (for example), if our operator is defined to be that differential operator restricted to Schwartz functions, then those solutions $u$ are in the kernel of the adjoint of the operator. In von Neumann's terms, this says that for $\Im(\lambda)>0$, the deficiency index is $2$. But for $\Im(\lambda)<0$ there are no $L^2$ solutions, so those deficiency indices are $0$. Thus, by von Neumann's criterion, there are no self-adjoint extensions of the original operator.
That is, we can define a symmetric operator $S$ by restricting the domain of the differential operator, but its adjoint $S^*$ will not by symmetric, and we cannot (by the previous paragraph) extend $S$ in any way whatsoever to achieve $S=S^*$. Thus, the adjoint will inevitably have some non-real eigenvalues.