I agree with the conclusion of user 'wondering', continuum states can indeed give a significant (often quite large) contribution to perturbation energies and wave functions. I would just like to give a few examples.
For the ground state of a perturbed hydrogen with Hamiltonian
$$
H=-\frac{1}{2}\nabla^2-\frac{1}{r}+W\equiv H^{(0)}+W \ ,
$$
and
$$
\psi^{(0)}_{100}(\vec{r})=2\exp(-r)Y_{00}(\theta,\phi) \ \ , \ \ E^{(0)}_1=-\frac{1}{2} \ ,
$$
the solution of the first-order Schrödinger equation
$$
(H^{(0)}-E^{(0)}_1)|\psi^{(1)}\rangle+(W-E^{(1)})|\psi^{(0)}_{100}\rangle=0
$$
can be written as an expansion over all (bound and scattering) hydrogenic eigenstates:
$$
\begin{aligned}
|\psi^{(1)}\rangle
&=
\frac{{\cal{P}}^\perp}{E^{(0)}_1-H^{(0)}}W|\psi^{(0)}_{100}\rangle \\
&=
\sum_{l=0}^{\infty}\sum_{m=-l}^{+l}
\left[
\sum_{\substack{n=l+1 \\ n\neq1}}^\infty\frac{\langle\psi^{(0)}_{nlm}|W|\psi^{(0)}_{100}\rangle}{E^{(0)}_1-E^{(0)}_n}|\psi^{(0)}_{nlm}\rangle+\int_0^\infty\mathrm{d}\varepsilon\frac{\langle\psi_{lm}^{(0)}(\varepsilon)|W|\psi^{(0)}_{100}\rangle}{E^{(0)}_1-\varepsilon}|\psi_{lm}^{(0)}(\varepsilon)\rangle
\right] \ ,
\end{aligned}
$$
where ${\cal{P}}^\perp=I-|\psi^{(0)}_{100}\rangle\langle\psi^{(0)}_{100}|$ is projecting onto the orthogonal complement space of $|\psi^{(0)}_{100}\rangle$, and intermediate normalization was assumed, so $\langle\psi^{(0)}_{100}|\psi^{(1)}\rangle=0$ and ${\cal{P}}^\perp|\psi^{(1)}\rangle=|\psi^{(1)}\rangle$. Zeroth-order states are orthonormal in the following sense:
$$
\begin{aligned}
\langle\psi^{(0)}_{nlm}|\psi^{(0)}_{n'l'm'}\rangle&=\delta_{nn'}\delta_{ll'}\delta_{mm'}
\ , \\
\langle\psi^{(0)}_{lm}(\varepsilon)|\psi^{(0)}_{l'm'}(\varepsilon')\rangle&=\delta(\varepsilon-\varepsilon')\delta_{ll'}\delta_{mm'}
\ , \\
\langle\psi^{(0)}_{nlm}|\psi^{(0)}_{l'm'}(\varepsilon')\rangle&=0 \ .
\end{aligned}
$$
The continuum part of the zeroth-order spectrum affects second- and third-order energies via the relations
$$
E^{(2)}=\langle\psi^{(0)}_{100}|W|\psi^{(1)}\rangle \ \ , \ \
E^{(3)}=\langle\psi^{(1)}|W|\psi^{(1)}\rangle-E^{(1)}\langle\psi^{(1)}|\psi^{(1)}\rangle \ .
$$
In the following, examples for $W$ are considered where $\psi^{(1)}$, $E^{(2)}$ and $E^{(3)}$ can be given in a closed analytical form, which makes it possible to assess the importance of continuum states in the sum-over-states expansions for each case. The bound state part of $E^{(2)}$ (denoted by $E^{(2)}_\text{b}$) will be calculated numerically, and compared to the analytically known $E^{(2)}=E^{(2)}_\text{b}+E^{(2)}_\text{c}$. We take $\lambda\geq0$ in all cases below.
$\bullet$ Case 1: $W={\cal{F}}z$
This is just the textbook case of the quadratic Stark effect in the ground state of hydrogen (see e.g. Ch. 33 of Schiff). From symmetry considerations, it is obvious that $E^{(1)}=0$ and that $\psi^{(1)}$ must have $p_z$ symmetry. After substituting the Ansatz $\psi^{(1)}(\vec{r})=\chi(r)Y_{10}(\theta,\phi)$, the first-order Schrödinger equation is reduced to an ordinary differential equation which is easily solved. We obtain
$$
\psi^{(1)}(\vec{r})=-\frac{2}{\sqrt{3}}{{\cal{F}}}r\left[1+\frac{r}{2}\right]\exp\left(-r\right)Y_{10}(\theta,\phi) \ ,
$$
from which the second-order energy is readily found:
$$
E=-\frac{1}{2}-\frac{9}{4}{\cal{F}}^2+{\cal{O}}({\cal{F}}^4) \ .
$$
The third-order energy is zero once again due to parity.
To find $E^{(2)}_\text{b}$, we need to compute the matrix elements of $z=r\sqrt{4\pi/3} \, Y_{10}(\theta,\phi)$:
$$
\begin{aligned}
\langle\psi^{(0)}_{n10}|z|\psi^{(0)}_{100}\rangle
&=
\frac{1}{\sqrt{3}}\int_0^\infty\mathrm{d}r \, r^2 \, R_{n1}(r) \, r \, R_{10}(r) \\
&=
\frac{1}{\sqrt{3n(n^2-1)}}\frac{192n^2}{(n+1)^5}\frac{(4)_{n-2}}{(n-2)!}\sum_{k=0}^{n-2}\frac{(-(n-2))_{k}(5)_{k}}{(4)_{k}k!}\left(\frac{2}{n+1}\right)^k \\
&=\frac{16}{\sqrt{3}}\left(\frac{n-1}{n+1}\right)^n\sqrt{\frac{n^7}{(n^2-1)^5}} \\
&\equiv f(n)\ ,
\end{aligned}
$$
where $(x)_n$ is the rising factorial. The above radial matrix element is a special case of a more general formula presented here. The sum over bound states is then calculated numerically:
$$
\frac{1}{{\cal{F}}^2}E^{(2)}_\text{b}=-\sum_{n=2}^\infty\frac{2n^2f^2(n)}{n^2-1}\approx-1.83163 \ .
$$
This implies $E^{(2)}_\text{c}/{\cal{F}}^2\approx-0.41837$, in good agreement with the value given in Ruffa, Am. J. Phys. 41 234 (1973); in other words, roughly $19\%$ of $E^{(2)}$ comes from the continuum. It may be tempting to explain this as a consequence of $H$ not having any exact bound states in the presence of a homogeneous electric field (only resonances); however, this cannot be true, as the same effect will be observed for Hamiltonians with genuine bound states.
$\bullet$ Case 2: $W=\lambda r$
Contrary to the previous case, the potential $-1/r+\lambda r$ (which is just a rescaled Cornell potential) gives rise purely to bound states. As far as I know, no closed form analytical solution can be found, but PT is easy to apply. We find $E^{(1)}=3\lambda/2$ and that $\psi^{(1)}$ must have $s$ symmetry; the first-order Schrödinger equation is again easily solved to yield
$$
\psi^{(1)}(\vec{r})=\lambda\left[3-r^2\right]\exp(-r)Y_{00}(\theta,\phi) \ ,
$$
$$
E=-\frac{1}{2}+\frac{3}{2}\lambda-\frac{3}{2}\lambda^2+\frac{27}{4}\lambda^3+{\cal{O}}(\lambda^4) \ .
$$
For e.g. $\lambda=0.01$, this perturbative result compares very favorably with the exact numerical ground state energy $E\approx-0.4851437$.
To compute the sum over bound states in $E^{(2)}$, we use another special case of the aforementioned general matrix element formula:
$$
\begin{aligned}
\langle\psi^{(0)}_{n00}|r|\psi^{(0)}_{100}\rangle
&=\frac{1}{\sqrt{n}}\frac{24n^2}{(n+1)^4}\frac{(2)_{n-1}}{(n-1)!}\sum_{k=0}^{n-1}\frac{(-(n-1))_{k}(4)_{k}}{(2)_{k}k!}\left(\frac{2}{n+1}\right)^k \\
&=\frac{3}{2}\delta_{n1}-8\left(\frac{n-1}{n+1}\right)^n\sqrt{\frac{n^5}{(n^2-1)^4}}(1-\delta_{n1}) \\
&\equiv f(n) \ .
\end{aligned}
$$
The sum is again computed numerically:
$$
\frac{1}{\lambda^2}E^{(2)}_\text{b}=-\sum_{n=2}^\infty\frac{2n^2f^2(n)}{n^2-1}\approx-1.074881 \ ,
$$
and it shows that roughly $28\%$ of the second-order energy comes from continuum states.
$\bullet$ Case 3: $W=\lambda/r^2$
The potential $-1/r+\lambda/r^2$ is just a rescaled Kratzer potential, for which the Schrödinger equation can be solved analytically (see e.g. Landau & Lifshitz III); the exact ground state energy reads
$$
\begin{aligned}
E
=-\frac{2}{\left(1+\sqrt{1+8\lambda}\right)^2}
=-\frac{1}{2}+2\lambda-10\lambda^2+56\lambda^3+{\cal{O}}(\lambda^4) \ .
\end{aligned}
$$
The PT series is convergent for $\lambda<1/8$. Of course, these corrections could also be found from the solution of the first-order Schrödinger equation:
$$
\psi^{(1)}(\vec{r})=4\lambda\left[r+\ln(2r)+\gamma-3\right]\exp(-r)Y_{00}(\theta,\phi) \ ,
$$
where $\gamma$ is the Euler-Mascheroni constant.
The bound part of $E^{(2)}$ is computed with the help of another special case of the aforementioned formula:
$$
\begin{aligned}
\langle\psi^{(0)}_{n00}|r^{-2}|\psi^{(0)}_{100}\rangle
&=\frac{1}{n^{3/2}}\frac{4}{n+1}\frac{(2)_{n-1}}{(n-1)!}\sum_{k=0}^{n-1}\frac{(-(n-1))_{k}}{(2)_{k}}\left(\frac{2}{n+1}\right)^k \\
&=-\frac{2}{n^{3/2}}\left[\left(\frac{n-1}{n+1}\right)^n-1\right]\\
&\equiv f(n) \ ,
\end{aligned}
$$
leading to the numerical sum
$$
\frac{1}{\lambda^2}E^{(2)}_\text{b}=-\sum_{n=2}^\infty\frac{2n^2f^2(n)}{n^2-1}\approx-1.56001 \ .
$$
This shows that roughly $84\%$ of the second-order energy is produced by the continuum.
$\bullet$ Case 4: $W=\lambda \delta(\vec{r})$
This is a confusing, ill-posed eigenvalue problem. Delta function interactions often appear as approximations to various processes in atomic physics (fine and hyperfine structure, vacuum polarization in bound state QED, finite nuclear size effects, etc.), to be used in first-order PT with the obvious result $E^{(1)}=\lambda/\pi$. However, the higher-order energy corrections all turn out to be divergent. The exact eigenvalues of this Hamiltonian are similarly divergent when attacked with Green's function techniques, but end up finite and "trivial" when calculated variationally, or in some regularization scheme: they simply reduce to the eigenvalues of the unperturbed, hydrogenic Hamiltonian! I wrote a very long (probably too long) answer about this anomaly and why one should not worry too much about it.
The divergence of $E^{(2)}$ and $E^{(3)}$ is immediately evident if one tries to compute them with
$$
\psi^{(1)}(\vec{r})=\frac{\lambda}{\pi}\left[2r+2\ln(2r)-\frac{1}{r}+2\gamma-5\right]\exp(-r)Y_{00}(\theta,\phi)
\ ;
$$
see the aforementioned answer for a derivation of $\psi^{(1)}$.
What makes this problem very deceptive is that all infinities of the PT come from the integral over continuum states, while the sum over bound states is finite:
$$
|\psi^{(0)}_{nlm}(0)|^2=\delta_{l0}\frac{1}{\pi n^3} \ \ , \ \
|\psi^{(0)}_{lm}(\varepsilon;0)|^2=\delta_{l0}\frac{1}{\pi}\frac{1}{\displaystyle1-\exp\left(-\frac{2\pi}{\sqrt{2\varepsilon}}\right)} \ ,
$$
$$
\begin{aligned}
\left(\frac{\pi}{\lambda}\right)^2E^{(2)}_\text{b}
&=-\sum_{n=2}^\infty\frac{2}{n(n^2-1)} \\
&=
-\sum_{n=2}^\infty\left[\frac{1}{n(n-1)}-\frac{1}{n(n+1)}\right] \\
&=-\frac{1}{2} \ ,
\end{aligned}
$$
$$
\begin{aligned}
\left(\frac{\pi}{\lambda}\right)^2E^{(2)}_\text{c}&=-\lim_{{\cal{K}}\rightarrow\infty}\int_0^{\cal{K}}\mathrm{d}\varepsilon\frac{2}{2\varepsilon+1}
\frac{1}{\displaystyle1-\exp\left(-\frac{2\pi}{\sqrt{2\varepsilon}}\right)} \\
&=-\lim_{{\cal{K}}\rightarrow\infty}\left[\frac{\sqrt{2{\cal{K}}}}{\pi}+\frac{1}{2}\ln(2{\cal{K}}+1)+\frac{1}{2}+{\cal{O}}({\cal{K}}^{-1/2})\right] \ .
\end{aligned}
$$
If one only considered the bound state part of the hydrogenic spectrum, then it would be easy to arrive at the wrong conclusion that $E^{(2)}$ is finite.
Continuum contributions to PT energies in the general case
The previous examples were special since the first-order wave function could be found analytically, and this took care of all zeroth-order state contributions. This is of course not possible in the general case, since usually we have only numerical solutions already for the zeroth-order problem. The cleanest way to obtain accurate PT energies is to reformulate the calculation as a variational problem: computing the exact second-order energy turns out to be equivalent to finding the minimum of the so-called Hylleraas functional:
$$
J[\phi]=\langle\phi|H^{(0)}-E^{(0)}|\phi\rangle+\langle\phi|W-E^{(1)}|\psi^{(0)}\rangle+\langle\psi^{(0)}|W-E^{(1)}|\phi\rangle \ .
$$
This optimization problem is not tied to sum-over-states expressions, and can be used to find highly accurate numerical PT energies. See Ch. 4.2 of Simple Theorems, Proofs, and Derivations in Quantum Chemistry by Mayer or Sec. II.a.25.$\beta$ of Quantum Mechanics of One- and Two-Electron Atoms by Bethe & Salpeter for further details.
One more example: the ground state energy of Helium
The prototypical application of the Hylleraas functional approach was the computation of the second-order energy of Helium and similar two-electron ions. The clamped nucleus Hamiltonian of such a system ($Z$-dependence scaled out with $\vec{r}{}'=Z\vec{r}$) reads
$$
\begin{aligned}
H
&=
-\frac{1}{2}\nabla_1^2-\frac{1}{2}\nabla_2^2-\frac{Z}{r_1}-\frac{Z}{r_2}+\frac{1}{|\vec{r}_1-\vec{r}_2|} \\
&=
Z^2\left[-\frac{1}{2}{\nabla'}_1^2-\frac{1}{2}{\nabla'}_2^2-\frac{1}{r'_1}-\frac{1}{r'_2}+\frac{1}{Z}\frac{1}{|\vec{r}{}'_1-\vec{r}{}'_2|}\right] \\
&\equiv
Z^2\left[H^{(0)}+W\right] \ ,
\end{aligned}
$$
the full electron-electron interaction being treated as a perturbation with "coupling constant" $1/Z$. This suggests that the nuclear charge-dependence of the energy is
$$
E(Z)=Z^2\sum_{k=0}^{\infty}\epsilon_k Z^{-k} \ ,
$$
and that $E^{(2)}$ is independent of $Z$. One can easily find $\epsilon_0$ and $\epsilon_1$, but the second- and higher-order energies can only be obtained numerically, and the above variational approach is required to get their accurate values:
$$
\epsilon_0=-1 \ \ , \ \ \epsilon_1=\frac{5}{8} \ \ , \ \
\epsilon_2\approx-0.157666429 \ \ , \ \ \epsilon_3\approx0.008699032 \ \ , \ \ ... \ .
$$
Using these coefficients, one can already approach the exact non-relativistic clamped-nucleus value $E\approx-2.903724$ for Helium.
I calculated the bound state part of $\epsilon_2$ similarly to the previous cases; without going into details, I found $\epsilon_{2,\text{b}}\approx-0.0631$, which is in good agreement with the old result of Scherr, J. Chem. Phys 33 317 (1960) and Knight, Scherr, Phys. Rev. 128 6 2675 (1962). This shows that only $40\%$ of the second-order energy is attributed to purely bound states, the other $60\%$ comes from states in which at least one of the electrons is in a scattering state.
