What does it mean for an operator to have both a continuous and discrete spectrum?

I'll try to give a mathematically rigorous answer which I hope will help you have a clearer understanding of the problem.

The setup

Let us consider the following hamiltonian in $L^2(\mathbb{R}^3)$ $$H = -\frac{\hbar^2}{2m}\nabla^2+V\tag{1}$$ and assume that $$V\in C^0(\mathbb{R}^3)\qquad |V(x)|< \frac{C}{|x|^a}\quad\text{with }C>0,\,a>2$$

Under this assumptions $H$ is self-adjoint on the domain of $H_0$, $D(H_0)$ and it's bounded from below. This is exactly the same problem you're facing (I choose it in $3D$, but the problem stays the same in $1$ and $2$ dim). The spectrum of $(1)$ is made up by a discrete part and an absolutely continuous part $$\sigma(H) = \sigma_p(H)\cup\sigma_{ac}(H) = \{E_1,\cdots, E_n\}\cup [0,+\inf)$$ Moreover we can decompose the Hilbert space as a direct product $L^2(\mathbb{R}^3) = \mathcal{H}_p\oplus\mathcal{H}_{ac}$.

The absence of positive eigenvalues means that the eigenvalue problem given by Shrödinger equation $$H\psi = E\psi \qquad E>0\tag{2}$$ has no solutions in $L^2(\mathbb{R}^3)$. But we can indeed construct a different class of solutions to $(2)$ called generalised eigenfunctions. This new class of solutions, together with the proper eigenfunctions, will help us understand the action of $H$ on a broader class of functions giving us the possibility to construct the discrete and continuous spectrum directly from them.

Building up the solution

In the case of $V=0$ we know that the generalised eigenfuncions are bounded solutions to the equation $$H_0\phi_0 = E\phi_0\qquad E>0$$ These solutions are given by plane-waves $$\phi_0(x,y) = e^{ikx}\qquad \forall k\in\mathbb{R}^3\quad k^2 = \frac{2mE}{\hbar^2}$$

Recall that from $\phi_0$ one can construct a unitary map $F$, which is the Fourier transform, that diagonalises the free hamiltonian operator $$[(FH_0F^{-1})\tilde{f}](k) = k^2 \tilde{f}(k)$$

Now we can generalise this procedure in the case of $V\neq 0$. So we have to study the above said equation $$H\phi= E\phi\qquad E>0$$ or, equivalently $$(-\nabla^2-k^2)\phi = -U\phi\qquad U = \frac{2mV}{\hbar^2}\tag{3}$$ To solve $(3)$ we look for a solution of the form $$\phi(x,t) = e^{ikx}+\eta(x,t)$$ where $\eta$ is a continuous, bounded function such that $\eta\to0$ for $|x|\to\infty$ which is what we would expect since whenever we go far enough away, we need to get back the free solution. Using the above form of $\phi$ in equation $(3)$ one gets an equation for $\eta$ $$(-\nabla^2-k^2)\eta = -U(e^{ik(\cdot)}+\eta)\tag{4}$$ We can write this solution in integral form since we know the green function for the free resolvent $(-\nabla^2-k^2)^{-1}$ using the limits $$\lim_{\epsilon\to 0}(-\nabla^2-k^2\pm i\epsilon)^{-1}$$ which are convenient since the operator $(-\nabla^2-k^2)$ is singular, and by so non invertible, whenever $k^2$ is an eigenvalue. Using this we find the integral form of $(4)$ $$\eta_\pm(x,k) = -\int\mathrm{d}y\,\frac{e^{\mp i|k||x-y|}}{4\pi|x-y|}U(y)\left(e^{iky}+\eta_\pm(y,k)\right)$$ If you want to know more about this specific result you can search for the integral kernel of the free resolvent.

Using this result we obtain $$\phi_\pm(x,k) = e^{ikx}-\int\mathrm{d}y\,\frac{e^{\mp i|k||x-y|}}{4\pi|x-y|}U(y)\phi_\pm(y,k)$$ this is called Lippman-Schwinger equation. What's important is that this equation $\forall k\in\mathbb{R}^3/\{0\}$ admits a unique solution such that $$\lim_{|x|\to\infty}\left(\phi_\pm(x,k)-e^{ikx}\right)=0$$ which is continuous and bounded and such that, for any $k$, is a solution to $$(-\nabla^2+U)\phi_\pm = k^2\phi_\pm$$

The action of the hamiltonian

Now that we have both generalised eigenfunctions and proper eigenfunctions we can diagonalise the hamiltonian. With the generalised eigenfuncions we can construct a generalised transform $F_\pm f = \tilde{f}_\pm\in L^2(\mathbb{R}^3)$ where $$(F_\pm f)(x) = \frac{1}{(2\pi)^{3/2}}\lim_{R\to\infty}\int_{|x|<R}\mathrm{d}x\,\overline{\phi_\pm(x,k)}f(x)\tag{5}$$ where the limit is there since we cannot define this transformation on $L^2(\mathbb{R}^3)$ without it, much like Fourier transform. It's important to note that $F_\pm :\mathcal{H}_{ac}\to L^2(\mathbb{R}^3)$ is unitary with inverse $F^*_\pm:L^2(\mathbb{R}^3)\to\mathcal{H}_{ac}$ given by $$(F^*_\pm f)(x) = \frac{1}{(2\pi)^{3/2}}\lim_{R\to\infty}\int_{|k|<R}\mathrm{d}k\,\phi_\pm(x,k)f(k)$$

At this point using result $(5)$ we can characterize the continuous spectrum of the hamiltonian $(1)$ using it's action on any state $f\in D(H)\cap \mathcal{H}_{ac}$ which is going to be $$(Hf)(x) = \frac{\hbar^2}{2m}\frac{1}{(2\pi)^{3/2}}\lim_{R\to\infty}\int_{|k|<R}\mathrm{d}k\,\phi_{\pm}(x,k)|k|^2\tilde{f}_\pm(k)$$

One can even study the evolution of a given state in the continuous part of the Hilbert space.

We can now glue together the generalised eigenfunctions, which give us a diagonalisation of the Hamiltonian in the absolutely continuous Hilbert space, with the proper eigenfunctions, which give us a diagonalisation of the hamiltonian in the discrete part of our Hilbert space, so to have a full diagonalisation of $H$ in the whole Hilbert space.

Left $\phi_n$ be a proper eigenfuncion of $H$, which is a solution to $H\phi_n = E_n\phi_n$ with $E_n<0$ and $\phi_n\in D(H)$. Moreover, denote $f_n = (\phi_n, f)$ the fourier coefficient of any $f\in L^2(\mathbb{R}^3)$. We consider the Hilbert space $L^2(\mathbb{R}^3)\oplus l^2$ whose generic element is a pair $\langle g, \{c_n\}\rangle$, with $g\in L^2(\mathbb{R}^3)$ and $\{c_n\}\in l^2$.

We define the following operator $$U_\pm : L^2(\mathbb{R}^3)\to L^2(\mathbb{R}^3)\oplus l^2\qquad U_\pm f = \langle \tilde{f}_\pm, \{f_n\}\rangle $$ which is a unitary operator such that $$(f,g) = \int\mathrm{d}k\,\overline{\tilde{f}_\pm(k)}\tilde{g}_\pm(k) + \sum_n f_n g_n$$ for any $f,g\in L^2(\mathbb{R}^3)$ and $\text{Ran}U_\pm = L^2(\mathbb{R}^3)\oplus l^2$. This helps us define a new scalar product on our whole Hilbert space and, by doing so, even the action of $H$ on the whole Hilbert space, which in fact will be given by $f\in D(H)$ $$(Hf)(x) = \frac{\hbar^2}{2m}\frac{1}{(2\pi)^{3/2}}\lim_{R\to\infty}\int_{|k|<R}\mathrm{d}k\,\phi_{\pm}(x,k)|k|^2\tilde{f}_\pm(k) +\sum_n E_n f_n\phi_n(x)$$ This gives a diagonalisation of $H$ over the whole Hilbert space, and by doing so it gives us the full spectrum, discrete and continuous. Moreover we can now precisely define the discrete and continuous Hilbert spaces as $$\mathcal{H}_{ac} = \{ f\in L^2(\mathbb{R}^3) |f_n = 0\;\forall n\}\qquad \mathcal{H}_{p} = \{ f\in L^2(\mathbb{R}^3) |\tilde{f}_\pm = 0\}$$

When it comes to operators acting on infinite dimensional Hilbert spaces, you need to specify some boundary condition before you find the eigenfunctions and eigenvalues. Playing around with any linear differential operator would show that different boundary conditions lead to different eigenvalues.

In the situation you're concerned with, the discrete spectrum is always associated with bound states, which are normalizeable. The continuous spectrum is not normalizeable. Normalizeability is a "boundary condition", so leads to different spaces of eigenvalues and eigenfunctions.

So there is nothing wrong with having both a discrete and continuous spectrum, as long as they are obtained by specifying different boundary conditions.

The spectrum of an operator acting on a separable Hilbert space (state space) can always be partitioned into three disjoint sets by the Lebesgue Decomposition Theorem. They are the pure point (p.p., i.e. eigenvalues), absolutely continuous (a.c.), and singular continuous (s.c.) parts of the spectrum. This wouldn’t come up in most linear algebra courses, but is very relevant to quantum mechanics. You gave an example of an operator with p.p. and a.c. parts. Singular continuous spectra are more exotic; the standard example is when the spectrum is a Cantor set. The only physical example I know of where the latter comes up (s.c.) is for the tight-binding (aka hopping) Hamiltonian associated with the Fibonacci quasicrystal; see Kohmoto, PRB 35,1020 (1987).

Most functional analysis texts would cover this. Here are some references. Under the title Functional Analysis, a classic, more mathematical book is by Yosida, a more application-heavy text is by Kreysig; then there’s the series of texts by Reed and Simon. A relevant operator-centered text is Perturbation Theory of Linear Operators by Kato. Granted, they each presume at least some background in undergrad real analysis.