What you do in perturbation theory is you assume the correct eigenvalue equation for a (hitherto) unknown correct wavefuction $|\psi_n\rangle$ and its associated eigenvalue $E_n$: $$ H|\psi_n\rangle = E_n |\psi_n\rangle,$$ where $H$ is the full Hamiltionian.
What you have is a main contribution $H_0$ (e.g. the Coulomb potential) and a perturbation $H'$ which is small compared to $H_0$ and that will therefore just slightly change the main solution.
You then assume that you can expand the correct solution for the energy and the wavefunction as a perturbative series, i.e. in terms that are smaller and smaller: $$ |\psi_n\rangle = |\psi_n\rangle^0 + |\psi_n\rangle^1 + |\psi_n\rangle^2 ...$$
$$ E_n = E_n^0 + E_n^1 + E_n^2 ... $$
where the exponents signify the order of the term. Higher order terms are smaller, and therefore only needed in you want higher precision.
Now, the full TISE becomes: $$ (H_0 + H')\,(|\psi_n\rangle^0 + |\psi_n\rangle^1 + |\psi_n\rangle^2 ... ) = (E_n^0 + E_n^1 + E_n^2 + ...)\,(|\psi_n\rangle^0 + |\psi_n\rangle^1 + |\psi_n\rangle^2+...)$$
Now you take the $0^{th}$ order equation -- where each term is of $O^{th}$ order:
$$H_0 |\psi_n\rangle^0\rangle = E_n^0 |\psi_n\rangle^0,$$
which is just the unperturbed equation, and therefore the starting point for any perturbative calculation.
Now look at the $1^{st}$ order equation -- remember that two $1^{st}$ order terms multiplied together give you a $2^{nd}$ order terms, whereas you only want to keep the $1^{st}$ order ones. $H'$ is first order:
$$ H_0 |\psi_n\rangle^1 + H'|\psi_n\rangle^0 = E_n^0|\psi_n\rangle^1 + E_n^1 |\psi_n\rangle^0.$$
What you are after is $E_n^1$, i.e. the $1^{st}$ order contribution to the energy.
Multiply by $^0\langle \psi_n|$ from the left:
$$ ^0\langle \psi_n|H_0 |\psi_n\rangle^1 + ^0\langle \psi_n|H'|\psi_n\rangle^0 = ^0\langle \psi_n|E_n^0|\psi_n\rangle^1 + ^0\langle \psi_n|E_n^1 |\psi_n\rangle^0,$$
$$E_n^0 \,^0\langle \psi_n|\psi_n\rangle^1 + \, ^0\langle \psi_n|H'|\psi_n\rangle^0 = E_n^0\,^0\langle \psi_n|\psi_n\rangle^1 + E_n^1 \,^0\langle \psi_n|\psi_n\rangle^0,$$
$$ \implies (E_n^0 - E_n^0) \,^0\langle \psi_n|\psi_n\rangle^1 + ^0\langle \psi_n|H'|\psi_n\rangle^0 = E_n^1 \,^0\langle \psi_n|\psi_n\rangle^0.$$
Assuming normalised states, $^0\langle \psi_n|\psi_n\rangle^0 = 1$, so:
$$ E_n^1 = ^0\langle \psi_n|H'|\psi_n\rangle^0 $$.
As other have noted, there's a mistake in you formula.
The same procedure applies for higher order corrections.
Best Answer
Here (see the file) is the discussion of the case you are interested in. It is referred as "nearly degenerate perturbation theory", and it is reduced to a usual degenerate perturbation theory by a simple rearrangement of terms in Hamiltonian. The resulting secular equation differs from the degenerate perturbation theory by a simple modification: \begin{equation} \begin{pmatrix} E_m^{(0)} + V_{mm} - \epsilon_\alpha & V_{mn}\\ V_{nm} & E_n^{(0)} + V_{nn} - \epsilon_\alpha \end{pmatrix} \begin{pmatrix} c_m^\alpha\\ c_n^\alpha \end{pmatrix} = 0. \end{equation} where $c^{\alpha}_k$ are the coefficients of the zero-order eigenstates $|\alpha\rangle = c^\alpha_m|\psi^{(0)}_m\rangle + c^\alpha_n|\psi_n^{(0)}\rangle$ with the first-order corrected energies $\epsilon_\alpha$.
The rearrangement is performed as follows. Consider a Hamiltonian $H = H_0 + V$, where \begin{equation} H_0 = \sum_i E_i^{(0)}|\psi_i^{(0)}\rangle\langle \psi_i^{(0)}|, \end{equation} \begin{equation} V = \sum_{ij} V_{ij}|\psi_i^{(0)}\rangle\langle \psi_j^{(0)}|, \end{equation} where $E_n^{(0)} \approx E_m^{(0)}$ for some $m$,$n$. Now, let me write \begin{equation} H = H_0' + V', \end{equation} where \begin{equation} H_0' = H_0 - \frac{E_m^{(0)} - E_n^{(0)}}{2}(|\psi_m^{(0)}\rangle\langle \psi_m^{(0)}| - |\psi_n^{(0)}\rangle\langle \psi_n^{(0)}|), \end{equation} \begin{equation} V' = V + \frac{E_m^{(0)} - E_n^{(0)}}{2}(|\psi_m^{(0)}\rangle\langle \psi_m^{(0)}| - |\psi_n^{(0)}\rangle\langle \psi_n^{(0)}|). \end{equation} The Hamiltonian $H'$ has two strictly degenerate eigenstates $|\psi_m^{(0)}\rangle$ and $|\psi_n^{(0)}\rangle$ with the energies $(E_m^{(0)} + E_n^{(0)})/2$, which allows to use usual degenerate perturbation theory by $V'$.