How can I compute variance of X^2 + XY where (X,Y) are normal distributed random variables? Given that $(X,Y)\sim N(0,0,\sigma_1^2,\sigma_2^2,\rho)$, I want to calculate $Var(X(X+Y))$. Here's what I've done so far:
$Var(X^2+XY)=E(X^2(Y+X)^2)-(E(X(Y+X))^2$
$=Cov(X^2,(X+Y)^2)+(Var(X)+(E(X))^2)(Var(X+Y)+(E(X+Y))^2)-[Cov(X,X+Y)+E(X)E(X+Y)]^2$
in which $Var(X),E(X)$ already known. $E(X+Y)$,$Var(X+Y),Cov(X,Y)$ are easy to know, but how can I simplify $Cov(X^2,(X+Y)^2)$?
Is there a simple method to calculate this?

By definition
\begin{align}
\operatorname{Var}(X(X+Y))
  &=E(X^2(X+Y)^2)-(E(X(X+Y)))^2
\\&=E(X^4)+2E(X^3Y)+E(X^2Y^2)-(E(X^2)+E(XY))^2. \tag{1}
\end{align}
Using Isserlis' theorem, in terms of the variances and covariance $\sigma_{11},\sigma_{22},\sigma_{12}$ of $X$ and $Y$,
\begin{align}
E(X^4)&=3\sigma_{11}^2
\\E(X^3Y)&=3\sigma_{11}\sigma_{12}    \tag{2}
\\E(X^2Y^2)&=\sigma_{11}\sigma_{22}+2\sigma_{12}^2. 
\end{align}
Substituting (2) into (1) and simplifying,
\begin{align}
\operatorname{Var}(X(X+Y))=2\sigma_{11}^2+4\sigma_{11}\sigma_{12}+\sigma_{11}\sigma_{22}+\sigma_{22}^2.
\end{align}
