Binomial Coefficient Recurrence Relation It turns out that,
$$
\sum_k \binom{m}{k}\binom{n}{k}\binom{m+n+k}{k} = \binom{m+n}{n} \binom{m+n}{m}
$$
where $\binom{m}{n}=0$ if $n>m$. One can run hundreds of computer simulations and this result always holds. Is there a mathematical proof for this?
Note 1: This sum is the discrete version of the following integral,
$$
\int_{-\infty}^{\infty}  \binom{m}{ax} \binom{n}{ax} \binom{m+n+ax}{ax}dx
=\frac{1}{a} \binom{m+n}{m}^2
$$
where $m,n \in \mathbb{N}$, $a \in (0,1)$ and,
$$
\binom{y}{x}:=\frac{\Gamma(1+y)}{\Gamma(1+x)\Gamma(1+y-x)}
$$
is the continuous version of the binomial coefficient.
This integral was posted here .
Note 2: Thanks to epi163sqrt for his brilliant application of  Egorychev's method. One question is: can this approach be extended to the continuous binomial coefficient? For the discrete binomial coefficient we have,
$$
\frac{1}{2\pi i}\oint_{\vert z \vert=1} \frac{(1+z)^k}{z^{j+1}} dz = \binom{k}{j}
$$
since,
$$
(1+z)^k = \sum_i \binom{k}{i} z^i
$$
and therefore $a_{-1}=\binom{k}{j}$. If one was to start with something like these two equations for the continuous binomial coefficient then Egorychev's method could be used for integrals.
 A: Note: Don Knuth et al. tell us in section 5.5 of Concrete Mathematics:
Binomial coefficients are like chameleons, changing their appearance easily.
Here are two variations based upon the same technique and underpinning the statement above. They might also be of interest for comparison. In fact these are special cases of instructive derivations answered by @MarkoRiedel some time ago. He has provided a great collection of worked out problems as pdf file (Egorychev) in his profile at this site.

We use the coefficient of operator $[z^k]$ to denote the coefficient of $z^k$ of a series. This way we can write for instance
\begin{align*}
\binom{m}{k}=[z^{k}](1+z)^m=[z^{m-k}]\frac{1}{(1-z)^{k+1}}\tag{1.1}
\end{align*}
where we also use the binomial series expansion and the binomial identity
\begin{align*}
\binom{-p}{q}=\binom{p+q-1}{q}(-1)^q\tag{1.2}
\end{align*}
Variation 1:

We obtain
\begin{align*}
\color{blue}{\sum_{k\geq 0}}&\color{blue}{\binom{m+n+k}{k}\binom{m}{k}\binom{n}{k}}\\
&=\sum_{k\geq 0}\binom{m}{k}[z^{n-k}](1+z)^n[w^{m+n}](1+w)^{m+n+k}\tag{2.1}\\
&=[z^n](1+z)^n[w^{m+n}](1+w)^{m+n}\sum_{k\geq 0}\binom{m}{k}\left(z(1+w)\right)^k\tag{2.2}\\
&=[z^n](1+z)^n[w^{m+n}](1+w)^{m+n}\left((1+z)+zw\right)^m\\
&=[z^n][w^{m+n}](1+w)^{m+n}\sum_{k=0}^m\binom{m}{k}(zw)^k(1+z)^{m+n-k}\\
&=\sum_{k=0}^m\binom{m}{k}[z^{n-k}](1+z)^{m+n-k}[w^{m+n-k}](1+w)^{m+n}\\
&=\sum_{k=0}^m\binom{m}{k}\binom{m+n-k}{n-k}\binom{m+n}{m+n-k}\tag{2.3}\\
&=\binom{m+n}{m}\sum_{k=0}^n\binom{m}{k}\binom{n}{n-k}\tag{2.4}\\
&\,\,\color{blue}{=\binom{m+n}{m}^2}
\end{align*}
and the claim follows.

Comment:

*

*In (2.1) we use $\binom{p}{q}=\binom{p}{p-q}$ and apply (1.1).


*In (2.2) we apply the rule $[z^{p-q}]A(z)=[z^p]z^qA(z)$ and do some rearrangements.


*In (2.3) we use the binomial identity
\begin{align*}
\binom{m+n-k}{n-k}\binom{m+n}{m+n-k}&=\frac{(m+n-k)!}{(n-k)!m!}\,\frac{(m+n)!}{(m+n-k)!k!}\\
&=\binom{m+n}{m}\binom{n}{n-k}
\end{align*}


*In (2.4) we apply Vandermonde's identity.
Variation 2:

We obtain
\begin{align*}
\color{blue}{\sum_{k\geq 0}}&\color{blue}{\binom{m+n+k}{k}\binom{m}{k}\binom{n}{k}}\\
&=\sum_{k\geq 0}\binom{m+n+k}{k}[z^{m-k}]\frac{1}{(1-z)^{k+1}}[w^{n-k}]\frac{1}{(1-w)^{k+1}}\tag{3.1}\\
&=[z^m][w^n]\frac{1}{(1-z)(1-w)}\sum_{k\geq 0}\binom{m+n+k}{k}\left(\frac{z}{1-z}\,\frac{w}{1-w}\right)^k\tag{3.2}\\
&=[z^m][w^n]\frac{1}{(1-z)(1-w)}\frac{1}{\left(1-\frac{zw}{(1-z)(1-w)}\right)^{m+n+1}}\tag{3.3}\\
&=[z^m][w^n](1-z)^{m+n}(1-w)^{m+n}\frac{1}{(1-z-w)^{m+n+1}}\\
&=[z^m]\frac{1}{1-z}[w^n](1-w)^{m+n}\frac{1}{\left(1-\frac{w}{1-z}\right)^{m+n+1}}\\
&=[z^m]\frac{1}{1-z}[w^n]\sum_{q\geq 0}\binom{m+n}{q}(-w)^q\\
&\qquad\cdot\sum_{r\geq 0}\binom{-(m+n+1)}{r}\left(-\frac{w}{1-z}\right)^r\tag{3.4}\\
&=[z^m]\frac{1}{1-z}[w^n]\sum_{q\geq 0}\binom{m+n}{q}(-w)^q\\
&\qquad\cdot\sum_{r\geq 0}\binom{m+n+r}{r}\left(\frac{w}{1-z}\right)^r\tag{3.5}\\
&=[z^m]\frac{1}{1-z}\sum_{j=0}^n\binom{m+n}{j}(-1)^j\binom{m+2n-j}{n-j}\left(\frac{1}{1-z}\right)^{n-j}\tag{3.6}\\
&=\sum_{j=0}^n\binom{m+n}{j}(-1)^j\binom{m+2n-j}{n-j}[z^m]\frac{1}{(1-z)^{n-j+1}}\\
&=\sum_{j=0}^n\binom{m+n}{j}(-1)^j\binom{m+2n-j}{n-j}\binom{-n+j-1}{m}(-1)^{m}\\
&=\sum_{j=0}^n\binom{m+n}{j}(-1)^j\binom{m+2n-j}{n-j}\binom{m+n-j}{m}\\
&=\binom{m+n}{m}\sum_{j=0}^n \binom{n}{j}(-1)^j\binom{m+2n-j}{n-j}\tag{3.7}\\
&=\binom{m+n}{m}\sum_{j=0}^n \binom{n}{j}(-1)^j[z^{n-j}](1+z)^{m+2n-j}\\
&=\binom{m+n}{m}[z^n](1+z)^{m+2n}\sum_{j=0}^n\binom{n}{j}(-1)^j\left(\frac{z}{1+z}\right)^j\\
&=\binom{m+n}{m}[z^n](1+z)^{m+2n}\left(1-\frac{z}{1+z}\right)^n\\
&=\binom{m+n}{m}[z^n](1+z)^{m+n}\\
&\,\,\color{blue}{=\binom{m+n}{m}^2}
\end{align*}
and the claim follows.

Comment:

*

*In (3.1) we use the coefficient of operator as in (1.1).


*In (3.2) we apply the rule $[z^{p-q}]A(z)=[z^p]z^qA(z)$.


*In (3.3) we apply 1.2 and use a geometric series expansion.


*In (3.4) we expand the series using the binomial series expansion.


*In (3.5) we use the binomial identity form (1.2).


*In (3.6) we select the coefficent of $w^n$.


*In (3.7) we use the binomial identity from (1.2).
