For this purpose ("beginning grad student") it would make sense to focus on the case that the discrete eigenvalues appear in an energy range that does not overlap with the continuous spectrum (say, $E>E_0$). One can then simply use the formulas for perturbation theory of discrete spectra, and replace the $\sum_n$ over energies $E>E_0$ in the continuous spectrum by an integral $\int_{E_0}^\infty dE\,\rho(E)$ weighted by the density of states $\rho(E)$.
To second order in $\lambda$ one thus has
$$E_n(\lambda) = E_n^{(0)} + \lambda \left \langle n^{(0)} \right |V\left |n^{(0)} \right \rangle + \lambda^2\sum_{k \ne n} \frac{\left |\left \langle k^{(0)} \right |V\left |n^{(0)} \right \rangle \right |^2} {E_n^{(0)} - E_k^{(0)}}$$
$$\qquad\qquad+ \lambda^2\int_{E_0}^\infty dE\,\rho(E) \frac{\left |\left \langle E \right |V\left |n^{(0)} \right \rangle \right |^2} {E_n^{(0)} - E}+{\cal O}(\lambda^3).$$
Here $|n^{(0)}\rangle$ and $|E\rangle$ are the unperturbed eigenstates in the discrete and continuous spectrum, respectively. I have assumed the eigenvalue $E_n^{(0)}$ has multiplicity one, otherwise the usual approach of degenerate perturbation theory applies.
Concerning the question, "how to approximate the continuous spectrum" in a numerical calculation, indeed one would discretize this. For example, if the continuous spectrum consists of plane waves $e^{ikx}$ one could impose periodic boundary conditions over a length $L$, and discretize $k_m=2\pi m/L$, $m\in\mathbb{Z}$. The integral $\int dE\,\rho(E)$ would then be a sum over $m$ and $E$ would be replaced by the dispersion relation $E(k_m)$.
The case of non-isolated discrete levels, embedded in the continuum, is more complicated. The coupling to the continuum turns the bound state into a resonance, a quasi-bound state with a finite life-time. This can be calculated in perturbation theory by a formula known as Fermi's golden rule. I have the impression from the way the question is formulated that this is not the problem at hand, correct me if I'm wrong.
