Binomial Sum Related to Fibonacci: $\sum\binom{n-i}j\binom{n-j}i=F_{2n+1}$ How would I prove
$$
\sum\limits_{\vphantom{\large A}i\,,\,j\ \geq\ 0}{n-i \choose j} {n-j \choose i}
=F_{2n+1}
$$
where $n$ is a nonnegative integer and $\{F_n\}_{n\ge 0}$ is a sequence of Fibonacci numbers?
Thank you very much! :)

Here is a solution using two complex variables.
Suppose we seek to evaluate in terms of Fibonacci numbers
$$\sum_{p,q\ge 0} {n-p\choose q} {n-q\choose p}.$$
We use the integrals
$${n-p\choose q} =
\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{(1-z)^{q+1} z^{n-p-q+1}} \; dz$$
and
$${n-q\choose p} = 
\frac{1}{2\pi i}
\int_{|w|=\epsilon} \frac{1}{(1-w)^{p+1} w^{n-p-q+1}} \; dw.$$
These correctly control the range
so we may let $p$ and $q$ go to infinity to get for the sum
$$\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{(1-z) z^{n+1}}
\frac{1}{2\pi i}
\int_{|w|=\epsilon} \frac{1}{(1-w) w^{n+1}}
\sum_{p,q\ge 0} \frac{z^{p+q} w^{p+q}}{(1-w)^p (1-z)^q}
\; dw \; dz
\\ = \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{(1-z) z^{n+1}}
\frac{1}{2\pi i}
\int_{|w|=\epsilon} \frac{1}{(1-w) w^{n+1}}
\\ \times \frac{1}{1-zw/(1-w)} \frac{1}{1-zw/(1-z)}
\; dw \; dz
\\ = \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n+1}}
\frac{1}{2\pi i}
\int_{|w|=\epsilon} \frac{1}{w^{n+1}}
\frac{1}{1-w-zw} \frac{1}{1-z-zw}
\; dw \; dz
\\ = \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n+2} (1+z)}
\frac{1}{2\pi i}
\int_{|w|=\epsilon} \frac{1}{w^{n+1}}
\frac{1}{w-1/(1+z)} \frac{1}{w-(1-z)/z}
\; dw \; dz.$$
We evaluate the inner integral using the fact that the residues of the
function in $w$ sum to zero. We  have two simple poles. We get for the
first pole at $w=(1-z)/z$ 
$$\frac{z^{n+1}}{(1-z)^{n+1}}
\frac{1}{(1-z)/z-1/(1+z)}
= \frac{z^{n+1}}{(1-z)^{n+1}}
\frac{z(1+z)}{(1-z)(1+z)-z}
\\ = \frac{z^{n+2}}{(1-z)^{n+1}}
\frac{1+z}{(1-z)(1+z)-z}.$$
Substituting this expression  into the outer integral we  see that the
pole at $z=0$ is canceled making for a contribution of zero.
For the second pole at $w=1/(1+z)$ we get
$$(1+z)^{n+1} \frac{1}{1/(1+z)-(1-z)/z}
= (1+z)^{n+1} \frac{z(1+z)}{z-(1-z)(1+z)}.$$
This yields the  contribution (taking into account the  sign flip from
the sum of residues)
$$\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n+2} (1+z)}
(1+z)^{n+1} \frac{z(1+z)}{1-z-z^2} \; dz
\\ = \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n+1}}
(1+z)^{n+1} \frac{1}{1-z-z^2} \; dz.$$
We  evaluate  this using  again  the fact  that  the  residues sum  to
zero. There are simple poles at $z=-\varphi$ and $z=1/\varphi.$
These yield
$$\left(\frac{1-\varphi}{-\varphi}\right)^{n+1} 
\frac{1}{-1+2\varphi} +
\left(\frac{1+1/\varphi}{1/\varphi}\right)^{n+1} 
\frac{1}{-1-2/\varphi}
\\= \frac{1}{\sqrt{5}} \frac{1}{\varphi^{2n+2}}
- \frac{1}{\sqrt{5}} \varphi^{2n+2}.$$
Taking into account the sign flip  this is obviously Binet / de Moivre
for $${\large F_{2n+2}}.$$
Remark. If we want to do this properly we also need to verify that
the residue  at infinity  in both  cases is zero.  For example  in the
first application we use the formula for the residue at infinity
$$\mathrm{Res}_{z=\infty} h(z)
= \mathrm{Res}_{z=0} 
\left[-\frac{1}{z^2} h\left(\frac{1}{z}\right)\right]$$
which in the present case gives for the inner term in $w$
$$- \mathrm{Res}_{w=0} \frac{1}{w^2} w^{n+1} 
\frac{1}{1/w-1/(1+z)} \frac{1}{1/w-(1-z)/z}
\\ = - \mathrm{Res}_{w=0}  w^{n+1} 
\frac{1}{1-w/(1+z)} \frac{1}{1-w(1-z)/z}$$
which is zero by inspection.
