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.