Dealing with a difficult sum of binomial coefficients, $\sum_{l=0}^{n}\binom{n}{l}^{2}\sum_{j=0}^{2l-n}\binom{l}{j} $ 
I am interested in finding an upper bound for the sum 
  $$F(n)= \sum_{l=0}^{n}\binom{n}{l}^{2}\;\sum_{j=0}^{2l-n}\binom{l}{j}$$ 

Ideally it should be possible to evaluate it exactly using some combinatorial identities or generating functions. Note that the $\binom{n}{l}$ term is squared, and that when $2l-n<0$, we define the inner to be zero.
Simple upper bound: By noticing that $2l-n\leq n$, we have that  $$\sum_{j=0}^{2l-n}\binom{l}{j}\leq 2^n,$$ and so $$F(n)\leq 2^n \sum_{l=0}^{n}\binom{n}{l}^{2}=2^n\binom{2n}{n}.$$ (See this question for a proof that $\sum_{l=0}^{n}\binom{n}{l}^{2}=\binom{2n}{n}$) Modifying the above based on when the inner sum is zero, that is when $2l\leq n$, we can replace this upper bound by $$F(n)\leq 2^{n-1}\binom{2n}{n}.$$ However it should be possible to do significantly better than this.
 A: Note: Here are some additional aspects which might be helpful.

If we replace in $F(n)$ the index $l$ with $n-l$ we obtain
  \begin{align*}
F(n)&=\sum_{l=0}^{n}{\binom{n}{l}}^2\,\sum_{j=0}^{2l-n}\binom{l}{j}\\
&=\sum_{l=0}^n\binom{n}{l}^2\,\sum_{j=0}^{n-2l}\binom{n-l}{j}\tag{1}\\
&=\sum_{l=0}^{\lfloor\frac{n}{2}\rfloor}\binom{n}{l}^2\,\sum_{j=0}^{n-2l}\binom{n-l}{n}
\end{align*}

Starting from (1) and considering the simple upper bound $n-l$ instead of $n-2l$ we obtain
\begin{align*}
\sum_{l=0}^n&\binom{n}{l}^2\,\sum_{j=0}^{n-2l}\binom{n-l}{j}\\
&\leq \sum_{l=0}^n\binom{n}{l}^2\,\sum_{j=0}^{n-l}\binom{n-l}{j}\\
&=\sum_{l=0}^n\binom{n}{l}^22^{l}\tag{2}\\
\end{align*}

Feeding OEIS with the first few numbers of (2)
  \begin{align*}
1,3,13,63,321,1683,8989
\end{align*}
  we get the sequence A001850 of Central Delannoy Numbers. These numbers have a generating function
  \begin{align*}
D(x)&=\frac{1}{\sqrt{1-6x+x^2}}\\
&=1+3x+13x^2+63x^3+321x^4+1683x^5+\cdots
\end{align*}
  from which an asymptotic formula for the coefficients $D_n$ can be derived (see p.2 in the reference above)
  \begin{align*}
D_n=\frac{(3+2\sqrt{2})^n}{\sqrt{\pi}\sqrt{\sqrt{2}-4}}\left(\frac{n^{-\frac{1}{2}}}{2}-\frac{23n^{-\frac{3}{2}}}{32(8+3\sqrt{2})}+\mathcal{O}(n^{-\frac{5}{2}})\right)
\end{align*}
Note: Although we get a nice generating function and an asymptotic formula for the Central Delannoy Numbers $D_n$ we can't find a closed expression for them in OEIS. This is a strong indicator, that there is none available. From this we may conclude that since OPs expression $F(n)$ is somewhat more complicated than $D_n$ no closed expression is available for $F(n)$ either.

Hints: 


*

*I suppose there is no appropriate generating function for $F(n)$ available, from which an asymptotic analysis of the coefficients similar to that of the $D_n$ could be done.

*Some more recherche regarding Delannoy numbers could provide additional insight.

*Maybe an asymptotic estimation using Stirlings formula together with some clever dissecting of the sums in $F(n)$ could be expedient.
A: By way  of an intermittent  progress report I submit  several integral
representations of the sum
$$\sum_{k=0}^n {n\choose k}^2 \times
\sum_{j=0}^{2k-n} {k\choose j}.$$
They were obtained using the Egorychev method.
The reader  is cordially invited  to comment, verify and  prove these.
Saddle point asymptotics  could be attempted if  these were univariate
integrals.   I present  them in  the  hope that  perhaps a  univariate
generating function can  be found where the catalog  of asymptotics of
coefficient  extractor integrals  could  be  applied. Sometimes  there
exist   substitutions    which   produce   radical    cancellation   /
simplification in these binomial sum integrals.

First representation:
$$\frac{1}{2\pi i}
\int_{|z|=\epsilon} 
\frac{(1+z)^{n}}{z^{n+1}} 
\frac{1}{2\pi i}
\int_{|w|=\epsilon} 
\frac{1}{1-w} \frac{1}{w^{n+1}}
(1 + w + w^2 z)^n
\; dw\; dz.$$

Second representation:
$$\frac{1}{2\pi i}
\int_{|z|=\epsilon} 
\frac{1}{z^{n+1}} (1+2z)^n (1+z)^n
\; dz
\\ + 
\frac{1}{2\pi i}
\int_{|z|=\epsilon} 
\frac{(1+z)^{n}}{z^{n+1}} 
\frac{1}{2\pi i}
\int_{|w|=\epsilon} 
\frac{1}{2w-1} \frac{1}{(1-w)^{n+1}}
\frac{1}{w^n}
((1-w)^2 z + w)^n
\; dw\; dz.$$

Third representation:
$$\frac{1}{2\pi i}
\int_{|z|=\epsilon} 
\frac{1}{z^{n}}
\frac{1}{2\pi i}
\int_{|w|=\epsilon} 
\frac{1}{w^{n+1}}\frac{(1+w)^n}{1-w}
\frac{1}{z-w/(1-w)}
((1+z)^2 w^2 + z)^n
\; dw\; dz.$$

Fourth representation:
$$\frac{1}{2\pi i}
\int_{|w|=\epsilon} 
\frac{1}{w^{n+1}}
\frac{1}{1-w}
\frac{1}{2\pi i}
\int_{|z|=\epsilon} 
\frac{1}{z^{n+1}}
(1 + z + wz)^n
(1 + w^2 z)^n
\; dz\; dw.$$

Fifth representation:
$$\frac{1}{2\pi i}
\int_{|w|=\epsilon} 
\frac{1}{w^{n+1}} \frac{1}{1-w}
\frac{1}{2\pi i}
\int_{|z|=\epsilon} 
\frac{1}{z^{n+1}}
\frac{1}{(1-z)^{n+1}}
(1 + w - wz)^n
(z + w^2 - w^2 z)^n
\; dz\; dw $$

Sixth representation:
$$\frac{1}{2\pi i}
\int_{|w|=\epsilon} 
\frac{1}{w^{n+1}} 
\frac{1}{1-w}
\frac{1}{2\pi i}
\int_{|z|=\epsilon} 
\frac{1}{z^{n+1}} \frac{1}{(1-z)^{n+1}}
(w+1-z)^n
\; dz\; dw 
\\ - \frac{1}{2\pi i}
\int_{|w|=\epsilon} 
\frac{1}{w^{2n}} 
\frac{1}{1-w}
\frac{1}{2\pi i}
\int_{|z|=\epsilon} 
\frac{1}{z^{n+1}} \frac{1}{(1-z)^{n+1}}
(w+1-z)^n (1-z + w^2 z)^n
\; dz\; dw .$$
A: First time contributor to Stack Exchange, so pardon my typesetting. My proof could use some more rigor, so I’ll only sketch the idea.
From an integral representation of the exponential generating function for the numbers F(n), I was able to do classic saddle point analysis.  I then used ‘de-Poissionization’ to extract an asymptotic formula for F(n).  I will give a closed form for only the most important entity, A, since the other factors are an awful combination of derivatives of a complicated function that is evaluated at the saddle point.
$$ F(n) \approx C \frac{A^{n+1}}{n+1} \bigl( 1 - \frac{A\, B}{n+2} \bigr) $$
$$ A=\bigl(\sqrt{y_0}+\sqrt{1+1/y_0}\bigr)^2 \, , y_0=\frac{1}{3}\bigl( -1 + \tau^{1/3} + \tau^{-1/3} \bigr) \, ,
\tau=\frac{1}{2}(25+3\sqrt{69}) .$$
Numerically,
$$ A=5.729031537980930838, B=4.1847320149263052, C=0.250939299689850514 .$$
As an example, F(2000)=1.03227E1513 and the approximation is F*(2000)=1.0318E1513, so that about 4 significant digits are obtained.  The proposer wanted a bound, and since the second term in the asymptotic series is negative, the first term is sufficient, for n large enough.
