Solution (source code)

= Solution

Let $H^{(0)}|n\rangle=E_n^{(0)}|n\rangle$ with a nondegenerate reference <eigenvalue>, and choose intermediate normalization $\langle n|n^{(1)}\rangle=0$. The first-order equation is
$$
(H^{(0)}-E_n^{(0)})|n^{(1)}\rangle
=-(H^{(1)}-E_n^{(1)})|n\rangle.
$$
Projection on $|j\rangle$, $j\ne n$, gives $\langle j|n^{(1)}\rangle=H_{jn}^{(1)}/(E_n^{(0)}-E_j^{(0)})$. The second-order projected equation gives $E_n^{(2)}=\langle n|H^{(1)}|n^{(1)}\rangle$, hence
$$
E_n^{(2)}
=-\sum_{j\ne n}
\frac{|H_{nj}^{(1)}|^2}{E_j^{(0)}-E_n^{(0)}}.
$$
The sum runs over a complete orthonormal reference basis except the reference state; continuum states require integration if present. Degenerate states with zero denominators require <degenerate perturbation theory> first. These are expansion coefficients in $E_n=E_n^{(0)}+\lambda E_n^{(1)}+\lambda^2E_n^{(2)}+\cdots$.