Since I had troubles with the notation the first time around, I did some real quantum chemistry (with a pencil and some paper) and finally was able to derive the whole RSPT ansatz. Since this is equivalent to the first approach, but yet also standalone, I decided to add it as a separate answer. In a few key steps, this approach is slightly different, so a one to one transformation will not be possible.
$$%Introducing some shortcuts
\require{cancel}
\newcommand{\op}[1]{\hat{\mathrm{#1}}} %\op{H}
\newcommand{\order}[1]{^{(#1)}} %E_n\order{1}
\newcommand{\overlap}[3]{\mathcal{#1}_{#2}\order{#3}} %\overlap{S}{m}{1}
\newcommand{\integral}[3]{\mathcal{#1}_{#2,#3}} %\integral{V}{i}{j}
\newcommand{\tagref}[1]{\mathrm{#1}}
$$
We start again by setting up the Hamiltonian, as a reference operator with the perturbation operator.
$$\op{H} = \op{H}\order{0} +\lambda\op{V} \tag1$$
The perturbed Schrödinger equation therefore becomes
$$(\op{H}\order{0} +\lambda\op{V})|n\rangle = E_n|n\rangle, \tag2$$
and we again develop the terms of the energy and wave function into Taylor power series.
\begin{align}
E_n &= E_n\order{0} + \lambda E_n\order{1} + \lambda^2 E_n\order{2} + \dots &
&= \sum_i \lambda^i E_n\order{i} \tag3\\
|n\rangle &= |n\order{0}\rangle + \lambda |n\order{1}\rangle
+ \lambda^2 |n\order{2}\rangle + \dots &
&= \sum_i \lambda^i |n\order{i}\rangle \tag4\\
\end{align}
We can write the unperturbed system, that provides us with a complete solution and therefore with a complete set of energies and wave functions, i.e. $n\in\mathbb{N}$.
$$\op{H}\order{0} |n\order{0}\rangle = E_n\order{0} |n\order{0}\rangle \tag5$$
We also know from this, that all solutions are orthonormal, i.e.
\begin{align}
\langle n\order{0} |n\order{0}\rangle &= 1,
\langle n\order{0} |m\order{0}\rangle &= 0. \tag6
\end{align}
Further we choose our perturbed wave function to be (intermediate) normalised, i.e.
\begin{align}
\langle n\order0|n\rangle &= 1, &
\langle n|n\rangle &= 1, &
\langle n\order{i}|n\order{j}\rangle &= \delta_{ij}
\begin{cases}
1, & i=j\\
0, & i\neq j\\
\end{cases}.\tag7
\end{align}
We continue with gathering the equations for all our orders.
\begin{align}
\lambda^0 &:&
\op{H}\order0 |n\order0\rangle
&= E_n\order0 |n\order0\rangle\tag{8a}\\
\lambda^1 &:&
\op{H}\order0 |n\order1\rangle + \op{V} |n\order0\rangle
&= E_n\order0 |n\order1\rangle + E_n\order1 |n\order0\rangle\tag{8b}\\
\lambda^2 &:&
\op{H}\order0 |n\order2\rangle + \op{V} |n\order1\rangle
&= E_n\order0 |n\order2\rangle + E_n\order1 |n\order1\rangle
E_n\order1 |n\order0\rangle\tag{8c}\\
\vdots\\
\lambda^n &:&
\op{H}\order0 |n\order{n}\rangle + \op{V} |n\order{n-1}\rangle
&= \sum_{i=0}^n E_n\order{i} |n\order{n-i}\rangle\tag{8d}\\
\end{align}
The two fundamental tricks of algebra are adding zero and multiplying by one.
We use this to set up our bases, we basically expand our bases into an auxiliary base (also known as 'resolution-of-the-identity'), which is possible, since the unperturbed system provides the complete set. This is what makes this approach somewhat different. While in my other answer we are just expanding the wafer function into an ordinary boring linear combination, we are using here the complete Hilbert space to do the same thing, but a little more elegant. However, this is also the key step: Multiplying by one.
At the same time we introduce an overlap integral $\integral{S}{m}{i}$, which makes writing a little easier further on. For any order $i\neq0$ we can therefore write:
\begin{align}
|n\order{i}\rangle &=
\sum_{m} |m\order0\rangle\langle m\order0|n\order{i}\rangle\\
&=
\sum_{m\neq n} |m\order0\rangle
\underbrace{\langle m\order0|n\order{i}\rangle}_{\overlap{S}{m}{i}}
+ |n\order0\rangle\cancelto{0}{\langle n\order0|n\order{i}\rangle}\\
|n\order{i}\rangle &= \sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{i}\\
\end{align}
The bases we need are for first and second order, therefore we write
\begin{align}
|n\order{1}\rangle &= \sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{1}, &
|n\order{2}\rangle &= \sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{2},
\end{align}
which seems rather trivial.
The next step is getting the expressions for the energies and corrections to the wave function.
For $\tagref{8a}$ this is trivial and it just yields the expectation value of the unperturbed Hamiltonian. We project on (multiply from the left with) $\langle n\order0|$, hence
\begin{align}
\langle n\order0|E_n\order0 |n\order0\rangle &=
\langle n\order0|\op{H}\order0 |n\order0\rangle \\
E_n\order0 \cancelto{1}{\langle n\order0|n\order0\rangle} &=
\langle n\order0|\op{H}\order0 |n\order0\rangle \\
E_n\order0 &= \langle n\order0|\op{H}\order0 |n\order0\rangle. \tag{9a}\\
\end{align}
For the first order correction to the energy we use the hermiticity of the Hamiltonian, i.e. we rewrite $\tagref{8a}$ in bra instead of the ket notation we used all way through.
$$\langle n\order0| \op{H}\order0 = E_n\order0 \langle n\order0|$$
We then further go about and project $\tagref{8b}$ on $\langle n\order0|$, hence
\begin{align}
\langle n\order0|\left(\op{H}\order0|n\order1\rangle
+ \op{V}|n\order0\rangle \right) &=
\langle n\order0|\left(E_n\order0|n\order1\rangle
+ E_n\order1|n\order0\rangle \right)\\
\langle n\order0|\op{H}\order0|n\order1\rangle
+ \langle n\order0| \op{V}|n\order0\rangle &=
\langle n\order0|E_n\order0|n\order1\rangle
+ \langle n\order0| E_n\order1|n\order0\rangle\\
E_n\order0 \cancelto{0}{\langle n\order0|n\order1\rangle}
+ \langle n\order0| \op{V}|n\order0\rangle &=
E_n\order0\cancelto{0}{\langle n\order0|n\order1\rangle}
+ E_n\order1 \cancelto{1}{\langle n\order0| n\order0\rangle}\\
E_n\order1 &= \langle n\order0| \op{V}|n\order0\rangle \tag{9b}
\end{align}
Now for the first order correction to the waffle function. I have to admit, this was the part where I struggled most. We will project on $\langle k\order0|\neq\langle n\order0|$, therfore obtaining:
\begin{align}
\langle k\order0|\left(\op{H}\order0|n\order1\rangle
+ \op{V}|n\order0\rangle \right) &=
\langle k\order0|\left(E_n\order0|n\order1\rangle
+ E_n\order1|n\order0\rangle \right)\\
\langle k\order0|\op{H}\order0|n\order1\rangle
+ \langle k\order0| \op{V}|n\order0\rangle &=
\langle k\order0|E_n\order0|n\order1\rangle
+ \langle k\order0| E_n\order1|n\order0\rangle\\
E_k\order0 \langle k\order0|n\order1\rangle
+ \langle k\order0| \op{V}|n\order0\rangle &=
E_n\order0\langle k\order0|n\order1\rangle
+ E_n\order1 \cancelto{0}{\langle k\order0| n\order0\rangle}\\
E_k\order0 \langle k\order0|\sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{1}
+ \langle k\order0| \op{V}|n\order0\rangle &=
E_n\order0\langle k\order0|\sum_{m\neq n}
|m\order0\rangle \overlap{S}{m}{1}\\
\end{align}
Before it gets to messy, we are going to reduce each of the terms individually and then reinsert, to obtain our expression.
\begin{align}
E_k\order0 \langle k\order0|\sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{1}
&= \sum_{m\neq n,k}E_k\order0
\cancelto{0}{\langle k\order0|m\order0\rangle} \overlap{S}{m}{1}
+ E_k\order0
\cancelto{1}{\langle k\order0|k\order0\rangle} \overlap{S}{k}{1}\\
&= E_k\order0 \overlap{S}{k}{1}\\
\end{align}
The second term stays the same, and the last one simplifies analogously.
\begin{align}
E_n\order0 \langle k\order0|\sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{1}
&= \sum_{m\neq n,k}E_n\order0
\cancelto{0}{\langle k\order0|m\order0\rangle} \overlap{S}{m}{1}
+ E_n\order0
\cancelto{1}{\langle k\order0|k\order0\rangle} \overlap{S}{k}{1}\\
&= E_n\order0 \integral{S}{k}{1}\\
\end{align}
And reinsterting we can write and resubstitute for the overlap integral $\overlap{S}{k}{1}$, with $m=k$:
\begin{align}
E_k\order0 \overlap{S}{k}{1} + \langle k\order0| \op{V}|n\order0\rangle
&= E_n\order0 \overlap{S}{k}{1}\\
\integral{S}{k}{1}
&= \frac{\langle k\order0| \op{V}|n\order0\rangle}{E_n\order0 -E_k\order0}\\
|n\order1\rangle
&= \sum_{k\neq n} |k\order0\rangle
\frac{\langle k\order0| \op{V}
|n\order0\rangle}{E_n\order0 -E_k\order0} \tag{10b}\\
\end{align}
To make matters a little bit easier further on, we will introduce another short notation for the perturbed integrals:
$$\integral{V}{i}{j} := \langle i\order0|\op{V}|j\order0\rangle$$
And therefore we can write
\begin{align}
|n\order1\rangle
&= \sum_{k\neq n} |k\order0\rangle
\frac{\integral{V}{k}{n}}{E_n\order0 -E_k\order0} \tag{10b}\\
\end{align}
For the second order correction to the energy we project $\tagref{8c}$ on $\langle n\order0|$. This is again fairly straight forward, the key is to insert $\tagref{10b}$ at the end to obtain the expression.
\begin{align}
\langle n\order0|\left(\op{H}\order0|n\order2\rangle
+ \op{V}|n\order1\rangle\right)
&= \langle n\order0|\left( E_n\order0|n\order2\rangle
+ E_n\order1|n\order1\rangle + E_n\order2|n\order0\rangle\right)\\
\langle n\order0|\op{H}\order0|n\order2\rangle
+ \langle n\order0|\op{V}|n\order1\rangle
&= E_n\order0 \langle n\order0|n\order2\rangle
+ E_n\order1 \langle n\order0|n\order1\rangle
+ E_n\order2 \langle n\order0|n\order0\rangle\\
\small
E_n\order0 \cancelto{0}{\langle n\order0|n\order2\rangle}
+ \langle n\order0|\op{V}|n\order1\rangle
&= \small
E_n\order0 \cancelto{0}{\langle n\order0|n\order2\rangle}
+ E_n\order1 \cancelto{0}{\langle n\order0|n\order1\rangle}
+ E_n\order2 \cancelto{1}{\langle n\order0|n\order0\rangle}\\
E_n\order2 &= \langle n\order0|\op{V}|n\order1\rangle\\
% &= \langle n\order0|\op{V} \sum_{k\neq n} |k\order0\rangle
% \frac{\langle k\order0| \op{V}
% |n\order0\rangle}{E_n\order0 -E_k\order0}\\
&= \sum_{k\neq n} \frac{\langle n\order0|\op{V} |k\order0\rangle
\langle k\order0| \op{V} |n\order0\rangle}{E_n\order0
-E_k\order0}\tag{9c}\\
&= \sum_{k\neq n} \frac{\integral{V}{n}{k}
\integral{V}{k}{n}}{E_n\order0 -E_k\order0}
\end{align}
Now let's get to the real fancy part: The second order correction to the basil. We project once again $\tagref{8c}$ on $\langle k\order| \neq \langle n\order 0|$. We do of course need our expression for the first order base and our expansion for the second order. But let us take one step at the time.
\begin{align}
\langle k\order0|\left( \op{H}\order0|n\order2\rangle
+ \op{V}|n\order1\rangle\right)
&= \langle k\order0 \left(
E_n\order0|n\order2\rangle
+ E_n\order1|n\order1\rangle
+ E_n\order2|n\order0\rangle\right)\\
\small
\langle k\order0| \op{H}\order0|n\order2\rangle
+ \langle k\order0| \op{V}|n\order1\rangle
&= \small
\langle k\order0| E_n\order0|n\order2\rangle
+ \langle k\order0| E_n\order1|n\order1\rangle
+ \langle k\order0| E_n\order2|n\order0\rangle\\
\end{align}
And we proceed again term by term, because it is pretty messy. We are starting with the right hand side, and there with the third term, as this one cancels, because of our normalisation.
\begin{align}
\langle k\order0| E_n\order2|n\order0\rangle
&= E_n\order2 \cancelto{0}{\langle k\order0| n\order0\rangle} =0\\
\end{align}
For the second term on the right we need our definition of the base. Since we are projecting on $k$ we choose the base to be represented by $m$. We also need the expression for the first order correction to the energy $\tagref{9b}$. Also note, that we are using the contracted notation.
\begin{align}
\langle k\order0| E_n\order1|n\order1\rangle
&= E_n\order1 \langle k\order0|n\order1\rangle\\
&= E_n\order1 \langle k\order0| \sum_{m\neq n} |m\order0\rangle
\frac{\integral{V}{m}{n}}{E_n\order0 - E_m\order0}\\
&= \sum_{m\neq n,k} E_n\order1 \cancelto{0}{\langle k\order0| m\order0\rangle}
\frac{\integral{V}{m}{n}}{E_n\order0 - E_m\order0}
+ E_n\order1 \cancelto{1}{\langle k\order0| k\order0\rangle}
\frac{\integral{V}{k}{n}}{E_n\order0 - E_k\order0}\\
&= \frac{\integral{V}{n}{n}\integral{V}{k}{n}}{E_n\order0 - E_k\order0}\\
\end{align}
For the first term on the right hand side we just need our expansion of the second order correction of the base.
\begin{align}
\langle k\order0| E_n\order0|n\order2\rangle
&= E_n\order0 \langle k\order0|n\order2\rangle\\
&= E_n\order0 \langle k\order0| \sum_{m\neq n} |m\order0\rangle
\overlap{S}{m}{2}\\
&= \sum_{m\neq n,k} E_n\order0 \cancelto{0}{\langle k\order0|m\order0\rangle}
\overlap{S}{m}{2}
+ E_n\order0 \cancelto{1}{\langle k\order0|k\order0\rangle}
\overlap{S}{k}{2}\\
&= E_n\order0 \overlap{S}{k}{2}\\
\end{align}
We have to use the same for the first term on the left hand side. In addition we are using the bra notation for the energy to simplify.
\begin{align}
\langle k\order0| \op{H}\order0|n\order2\rangle
&= E_k\order0 \langle k\order0|n\order2\rangle\\
&= E_k\order0 \langle k\order0|\sum_{m\neq n} |m\order0\rangle
\overlap{S}{m}{2} \\
&= \sum_{m\neq n,k} E_k\order0 \cancelto{0}{\langle k\order0|m\order0\rangle}
\overlap{S}{m}{2}
+ E_k\order0 \cancelto{1}{\langle k\order0|k\order0\rangle}
\overlap{S}{k}{2}\\
&= E_k\order0 \overlap{S}{k}{2}\\
\end{align}
For the second term on the left side we once again need the expression for the first order correction to bass, choosing $m$, as we are already projecting on $k$.
\begin{align}
\langle k\order0| \op{V}|n\order1\rangle
&= \langle k\order0| \op{V} \sum_{m\neq n}|m\order0\rangle
\frac{\integral{V}{m}{n}}{E_n\order0 -E_m\order0}\\
&= \sum_{m\neq n}\frac{\integral{V}{k}{m}
\integral{V}{m}{n}}{E_n\order0 -E_m\order0}\\
\end{align}
Now, let's put all the stuffing back into the turkey. Then we rearrange to complete bird and solve for our desired expression for the base.
\begin{align}
E_k\order0 \overlap{S}{k}{2}
+ \sum_{m\neq n}\frac{\integral{V}{k}{m}
\integral{V}{m}{n}}{E_n\order0 -E_m\order0}
&= E_n\order0 \overlap{S}{k}{2}
+ \frac{\integral{V}{n}{n}\integral{V}{k}{n}}{E_n\order0 - E_k\order0}\\
(E_n\order0 - E_k\order0)\overlap{S}{k}{2}
&= \sum_{m\neq n}\frac{\integral{V}{k}{m}
\integral{V}{m}{n}}{E_n\order0 -E_m\order0}
- \frac{\integral{V}{n}{n}\integral{V}{k}{n}}{E_n\order0 - E_k\order0}\\
\overlap{S}{k}{2}
&= \sum_{m\neq n}\frac{\integral{V}{k}{m}
\integral{V}{m}{n}}{(E_n\order0 -E_m\order0)(E_n\order0 - E_k\order0)}
- \frac{\integral{V}{n}{n}\integral{V}{k}{n}}{(E_n\order0 - E_k\order0)^2}\\
|n\order2\rangle
&= \sum_{k\neq n} \sum_{m\neq n} |k\order0\rangle
\frac{\integral{V}{k}{m}
\integral{V}{m}{n}}{(E_n\order0 -E_m\order0)(E_n\order0 - E_k\order0)}\\
%split into two lines
&\phantom{=~}
- \sum_{k\neq n} |k\order0\rangle
\frac{\integral{V}{k}{n}\integral{V}{n}{n}}{(E_n\order0 - E_k\order0)^2}\\
\end{align}
And last, but not least, let's remove that nasty shorthand and make everything pretty again.
\begin{align}
|n\order2\rangle
&= \sum_{k\neq n} \sum_{m\neq n} |k\order0\rangle
\frac{
\langle k\order0|\op{V}|m\order0\rangle
\langle m\order0|\op{V}|n\order0\rangle
}{(E_n\order0 -E_m\order0)(E_n\order0 - E_k\order0)}
- \sum_{k\neq n} |k\order0\rangle
\frac{
\langle k\order0|\op{V}|n\order0\rangle
\langle n\order0|\op{V}|n\order0\rangle
}{(E_n\order0 - E_k\order0)^2}\\
\end{align}
And there we are. Pretty.
Now you will probably have noticed, that you were one step away from the ultron goal. In your last equation you simply had to reduce to the non-zero parts.
\begin{equation}
( E_n^0 - E_k^0 ) \langle k^0 |n^2\rangle = \sum_{m \neq n} \frac{ \langle k^0| \hat{V} |m^0 \rangle \langle m^0 | \hat{V} | n^0 \rangle}{E_n^0 - E_m^0} - \langle n^0 | \hat{V} | n^0 \rangle \sum_{m \neq n} \frac{\langle k^0 |m^0 \rangle \langle m^0 | \hat{V} | n^0 \rangle}{E_n^0 - E_m^0}
\end{equation}
\begin{align}
( E_n^0 - E_k^0 ) \langle k^0 |n^2\rangle
&= \sum_{m \neq n}
\frac{
\langle k^0| \hat{V} |m^0 \rangle
\langle m^0| \hat{V} |n^0 \rangle
}{E_n^0 - E_m^0}\\
&\phantom{=~}
- \sum_{m \neq n,k}
\frac{ \cancelto{0}{\langle k^0 |m^0 \rangle}
\langle n^0 | \hat{V} | n^0 \rangle
\langle m^0 | \hat{V} | n^0 \rangle
}{E_n^0 - E_m^0}\\
&\phantom{=~}
- \frac{ \cancelto{1}{\langle k^0 |k^0 \rangle}
\langle k^0 | \hat{V} | n^0 \rangle
\langle n^0 | \hat{V} | n^0 \rangle
}{E_n^0 - E_k^0}\\
( E_n^0 - E_k^0 ) \langle k^0 |n^2\rangle
&= \sum_{m \neq n}
\frac{
\langle k^0| \hat{V} |m^0 \rangle
\langle m^0| \hat{V} |n^0 \rangle
}{E_n^0 - E_m^0}
- \frac{
\langle k^0 | \hat{V} | n^0 \rangle
\langle n^0 | \hat{V} | n^0 \rangle
}{E_n^0 - E_k^0}\\
\langle k^0 |n^2\rangle
&= \sum_{m \neq n}
\frac{
\langle k^0| \hat{V} |m^0 \rangle
\langle m^0| \hat{V} |n^0 \rangle
}{(E_n^0 - E_m^0)(E_n^0 - E_k^0)}
- \frac{
\langle k^0 | \hat{V} | n^0 \rangle
\langle n^0 | \hat{V} | n^0 \rangle
}{\left(E_n^0 - E_k^0\right)^2}\\
\end{align}
Now you just have to expand into the complete base by multiplying with $\sum_{k\neq n}|k^0\rangle$ and you get to the same expression as above.
And as a theoretician a proof is only complete with a coffee stain, and trust me, there have been many.
Image is courtesy of Roger Karlsson (http://www.free-photo-gallery.org/photos/coffee-stain/) obtained from flicker.
