If we approach this problem as you state it, without assuming any context, we have the following: we are given a matrix in an unspecified basis, which Cosmas identified as the energy eigenstates of a harmonic oscillator, and which we are told corresponds to the position operator.
In this (Hilbert) basis, the delta functions are the standard basis elements, which are not position eigenfunctions, since they are not eigenfunctions of the given matrix.
Let's call our basis elements $|0\rangle,|1\rangle,|2\rangle,\ldots,|k\rangle,\ldots$. Their linear span is dense in the state space, and a general element can be written as a limit w.r.t. the obvious inner product that makes this into an orthonormal basis, which is a square integrable sequence $(a_k)_{k\ge0} = (a_0,a_2,a_3,\ldots)$.
Your question is how to find the eigenvalues if you cannot write down a characteristic polynomial. When the matrix has a nice form like the present one, in which the matrix entries are an expression in number $n$ of the row they're in, we get a recurrence that allows you to find the coefficients one by one.
In our example, for $(a_k)$ to form an eigenvector with eigenvalue $\lambda$ we must have
$$\sqrt1a_1 = \lambda a_0$$
and
$$\sqrt{n}a_{n-1} + \sqrt{n+1}a_{n+1} = \lambda a_n.$$
This defines a recurrence that can be simplified by writing $b_k = \sqrt{k!}a_k$, for which we get
$$b_1 = \lambda b_0$$
and
$$nb_{n-1} + b_{n+1} = \lambda b_n.$$
The solution with $b_0 = 0$ gives the $0$ vector, so we can assume that $b_0 = 1$. It is clear that for every $\lambda$ we get a sequence $b_k$ that is a solution ($\lambda = 0$ gives Cosmas' solution, as it should). This gives an infinite series in the original basis, whose limit, if it existed, would be the position eigenstate for the position $\lambda$. Note that going back to the harmonic oscillator, this is an explicit sum in energy eigenstates whose explicit form is known.
What is not so clear is if the associated sequence $(a_k)$ is actually square integrable, i.e. if $\langle a,a\rangle = \sum a_k^2 = \sum \frac{b_k^2}{k!} < \infty$, for which we need to bound the growth of the $b_k$. Even better would be to find a closed form expression for the $b_k$.
We can start by setting $f(x) = \sum_{k = 0}^\infty b_kx^k$ be the generating function of the $b_k$. Then the recurrence relation implies that $f - \lambda xf + x^2(xf)' = f - \lambda xf + x^2f + x^3f'$ is a power series whose coefficient of $x^{n+1}$ is $nb_{n-1} + b_{n+1} - \lambda b_n$ for $n\ge 0$, so this is 0, while the constant term is $b_0 = 1$, so finally we can obtain
$$(1 - \lambda x + x^2)f(x) + x^3f'(x) = 1.$$
This is general as well. In our case, this doesn't look very easy at all, but sympy actually finds the solution
$$f{\left (x \right )} = \frac{1}{x} \left(C_{1} + \int \frac{e^{\frac{\lambda}{x}}}{x^{2}} e^{- \frac{1}{2 x^{2}}}\, dx\right) e^{- \frac{1}{2 x^{2}} \left(2 \lambda x - 1\right)}$$
For nice $f$, you might find a closed form expression for the $b_n$, hence the $a_n$ from its Taylor coefficients. In this case, if anyone knows how to find a closed form expression for the coefficients from this (or in any other way) or how to bound the coefficients $b_n$, please comment.
As a final remark, the sum of the series may exist in some formal sense even if it lies outside the state space. If the sum of the squares diverges, that means that it isn't normalizable, so that the standard probabilistic interpretation of the state vector will present some difficulties. As ACuriousMind remarked, the fact that there is an uncountable number of position eigenvalues should imply that none of these sums actually is square integrable.