All sums below are by non-negative $i,j$ unless opposite is stated explicitly.
Denote $h(x,y)=\sum \binom{i+j}i x^iy^j=\frac1{1-(x+y)}$, $f(x,y)=\sum \binom{i+j}i^2 x^iy^j$. We want to prove that $2xyf^2(x^2,y^2)$ is an odd (in both $x$ and $y$) part of the function $h(x,y)$. In other words, we want to prove that $$2xyf^2(x^2,y^2)=\frac14\left(h(x,y)+h(-x,-y)-h(x,-y)-h(-x,y)\right)=\frac{2xy}{1-2(x^2+y^2)+(x^2-y^2)^2}.$$ So, our identity rewrites as $$f(x,y)=(1-2(x+y)+(x-y)^2)^{-1/2}=:f_0(x,y)$$ This is true for $x=0$, both parts become equal to $1/(1-y)$. Next, we find a differential equation in $x$ satisfied by the function $f_0$. It is not a big deal: $$\left(f_0(1-2(x+y)+(x-y)^2)\right)'_x=(x-y-1)f_0.$$ Since $f_0(0,y)$ and this equation uniquely determine the function $f_0$, it remains to check that this holds for $f(x,y)$, which is a straightforward identity with several binomials.
(Warning: I did not check it myself, but if your statement, and everything written by me above, is correct, this identity must be true.)