I deem it would be of general interest
to develop this scheme in a step-by-step fashion for the more general case in which
$k$ extractions are made from an urn containing $n$ red molecules at start ($1 \leqslant n$),
and inserting one blue molecule after each extraction.
We can then construct the $\left( {0,\; \ldots ,\;n} \right) \times \left( {0,\; \ldots ,\;n} \right)$ transition matrix, giving the probability of getting
$m(k)$ blue molecules at the $k$-th extraction, given that there are $m(k-1)$
blue molecules resulting from precedent extractions, which will clearly be:
$$
\overline {\mathbf{T}_{\,\mathbf{n}} } = \;\left[ {\begin{array}{*{20}c}
{m_{\,k - 1} \backslash m_{\,k} } &| & 0 & 1 & 2 & \cdots & {n - 1} & n \\
\hline
0 &| & {0/n} & {n/n} & {} & {} & {} & {} \\
1 &| & {} & {1/n} & {\left( {n - 1} \right)/n} & {} & {} & {} \\
2 &| & {} & {} & {2/n} & \ddots & {} & {} \\
\vdots &| & {} & {} & {} & \ddots & {} & {} \\
{n - 1} &| & {} & {} & {} & {} & {\left( {n - 1} \right)/n} & {1/n} \\
n &| & {} & {} & {} & {} & {} & {n/n} \\
\end{array} } \right]
$$
where empty cells contain $0$, and
$\mathbf{T}_{\,\mathbf{n}} = \overline {\overline {\mathbf{T}_{\,\mathbf{n}} } } = \text{transpose}\left( {\overline {\mathbf{T}_{\,\mathbf{n}} } } \right)$
So, representing the probabilities of having $0,\; \ldots ,\;n$ blue molecules by the column vector $\mathbf{m}_{\,\mathbf{n}} \left( k \right)$, we can write:
$$
\mathbf{m}_{\,\mathbf{n}} \left( k \right) = \mathbf{T}_{\,\mathbf{n}} \;\mathbf{m}_{\,\mathbf{n}} \left( {k - 1} \right) = \mathbf{T}_{\,\mathbf{n}} ^{\,\mathbf{k}} \;\mathbf{m}_{\,\mathbf{n}} \left( 0 \right) = \mathbf{T}_{\,\mathbf{n}} ^{\,\mathbf{k}} \;\left[ {\begin{array}{*{20}c}
1 \\
0 \\
\vdots \\
0 \\
\end{array} } \right]
$$
The interesting fact is that the matrix $\mathbf{T}_{\,\mathbf{n}} $ diagonalizes nicely into:
$$
\mathbf{T}_{\,\mathbf{n}} = \mathbf{C}_{\,\mathbf{n}} \;\left( {\left( {i/n} \right) \circ \mathbf{I}_{\,\mathbf{n}} } \right)\;\mathbf{C}_{\,\mathbf{n}} ^{\, - \,\mathbf{1}} = \mathbf{C}_{\,\mathbf{n}} \;\left( {\left( {i/n} \right) \circ \mathbf{I}_{\,\mathbf{n}} } \right)\;\mathbf{C}_{\,\mathbf{n}}
$$
where the matrix in the centre indicates the diagonal matrix
$$
\left( {\left( {i/n} \right) \circ \mathbf{I}_{\,\mathbf{n}} } \right) = \left[ {i/n\;\delta _{\,i,\,j} } \right]_{\,\mathbf{n}}
$$
while the matrix $\mathbf{C}_{\,\mathbf{n}} $ represents
$$
\mathbf{C}_{\,\mathbf{n}} = \mathbf{C}_{\,\mathbf{n}} ^{\, - \,\mathbf{1}} = \left[ {\left( { - 1} \right)^{\,i} \left( \begin{gathered}
n - j \\
n - i \\
\end{gathered} \right)} \right]_{\,\mathbf{n}}
$$
Therefore
$$
\begin{gathered}
\mathbf{m}_{\,\mathbf{n}} \left( k \right) = \mathbf{T}_{\,\mathbf{n}} ^{\,\mathbf{k}} \;\mathbf{m}_{\,\mathbf{n}} \left( 0 \right) = \mathbf{C}_{\,\mathbf{n}} \;\left( {\left( {i/n} \right) \circ \mathbf{I}_{\,\mathbf{n}} } \right)^{\,\mathbf{k}} \;\mathbf{C}_{\,\mathbf{n}} ^{\, - \,\mathbf{1}} \;\mathbf{m}_{\,\mathbf{n}} \left( 0 \right) = \hfill \\
= \mathbf{C}_{\,\mathbf{n}} \;\left( {\left( {i/n} \right)^{\,k} \circ \mathbf{I}_{\,\mathbf{n}} } \right)\;\mathbf{C}_{\,\mathbf{n}} \;\mathbf{m}_{\,\mathbf{n}} \left( 0 \right) = \hfill \\
= \left[ {\sum\limits_{0\, \leqslant \,l\, \leqslant \,n} {\left( { - 1} \right)^{\,i} \left( \begin{gathered}
n - l \\
n - i \\
\end{gathered} \right)\left( {l/n} \right)^{\,k} \left( { - 1} \right)^{\,l} \left( \begin{gathered}
n \\
n - l \\
\end{gathered} \right)} } \right]_{\,\mathbf{n}} = \hfill \\
= \left[ {\sum\limits_{0\, \leqslant \,l\,\left( { \leqslant \,i} \right)} {\left( { - 1} \right)^{\,i - l} \left( \begin{gathered}
n \\
n - i \\
\end{gathered} \right)\left( \begin{gathered}
i \\
l \\
\end{gathered} \right)\left( {l/n} \right)^{\,k} } } \right]_{\,\mathbf{n}} = \hfill \\
= \left[ {\frac{1}
{{n^{\,k} }}\left( \begin{gathered}
n \\
i \\
\end{gathered} \right)\sum\limits_{0\, \leqslant \,l\,\left( { \leqslant \,i} \right)} {\left( { - 1} \right)^{\,i - l} \left( \begin{gathered}
i \\
l \\
\end{gathered} \right)l^{\,k} } } \right]_{\,\mathbf{n}} = \hfill \\
= \left[ {\frac{{i!}}
{{n^{\,k} }}\left( \begin{gathered}
n \\
i \\
\end{gathered} \right)\left\{ \begin{gathered}
k \\
i \\
\end{gathered} \right\}} \right]_{\,\mathbf{n}} = \left[ {\frac{{n!}}
{{n^{\,k} \left( {n - i} \right)!}}\left\{ \begin{gathered}
k \\
i \\
\end{gathered} \right\}} \right]_{\,\mathbf{n}} \hfill \\
\end{gathered}
$$
Now, if after $k$ extractions with replacement, ending with $m_k$ blue molecules,
we extract (without replacement) all the $n$ balls from the urn and ask for the probability of the last to be red, that will be:
$$
\frac{{p(k,n)}}
{{p\left( {m_{\,k} } \right)}} = \frac{{\left( \begin{gathered}
n - 1 \\
m_{\,k} \\
\end{gathered} \right)}}
{{\left( \begin{gathered}
n \\
m_{\,k} \\
\end{gathered} \right)}} = \frac{{n - m_{\,k} }}
{n}\quad \left| \begin{gathered}
\;0 \leqslant m_{\,k} \leqslant n \hfill \\
\;1 \leqslant n \hfill \\
\end{gathered} \right.
$$
So, finally we arrive at:
$$
\begin{gathered}
p(k,n) = \frac{1}
{{n^{\,k} }}\sum\limits_{0\, \leqslant \,l\,\left( { \leqslant \,n} \right)} {\left( {\frac{{n - l}}
{n}} \right)\frac{{n!}}
{{\left( {n - l} \right)!}}\left\{ \begin{gathered}
k \\
l \\
\end{gathered} \right\}} = \hfill \\
= \frac{1}
{{n^{\,k} }}\sum\limits_{0\, \leqslant \,l\,\left( { \leqslant \,n} \right)} {\left\{ \begin{gathered}
k \\
l \\
\end{gathered} \right\}\left( {n - 1} \right)^{\,\underline {\,l\,} } } = \left( {\frac{{n - 1}}
{n}} \right)^{\,k} \hfill \\
\end{gathered}
$$
which confirms and generalize the formula already given.
