Solution (source code)

= Solution

Project the first-order equation onto $\Psi_j^{(0)}$, with $j\ne i$. By <orthonormality> and the <Hermitian operator> property,
$$
\begin{aligned}
\langle\Psi_j^{(0)}|(\hat H^{(0)}-E_i^{(0)})|\Psi_i^{(1)}\rangle
&=(E_j^{(0)}-E_i^{(0)})c_{ij},\\
\langle\Psi_j^{(0)}|(\hat H^{(1)}-E_i^{(1)})|\Psi_i^{(0)}\rangle
&=H_{ji}.
\end{aligned}
$$
Consequently
$$
\boxed{c_{ij}=\frac{H_{ji}}{E_i^{(0)}-E_j^{(0)}}
=\frac{H_{ij}^*}{E_i^{(0)}-E_j^{(0)}}\qquad(j\ne i),}
$$
and
$$
\boxed{\Psi_i^{(1)}=\sum_{j\ne i}
\frac{\langle\Psi_j^{(0)}|\hat H^{(1)}|\Psi_i^{(0)}\rangle}
{E_i^{(0)}-E_j^{(0)}}\Psi_j^{(0)}.}
$$
The numerator order is important: the question defines $H_{ij}$ with $i$ on the bra, whereas projecting the correction onto state $j$ gives $H_{ji}$. The two coincide only when the relevant <matrix element> is real.

Excluding an $i$ component chooses $\langle\Psi_i^{(0)}|\Psi_i^{(1)}\rangle=0$, a permitted first-order phase and <wavefunction normalization> convention. A component parallel to $\Psi_i^{(0)}$ is not determined by the projected equation. Nondegeneracy ensures nonzero denominators; a small <energy> gap can still make the perturbative expansion inaccurate unless $|\lambda H_{ji}|\ll|E_i^{(0)}-E_j^{(0)}|$.