# Covariance of Bivariate Gaussian RV which is a function of 2 other Gaussian RVs

If $X \sim N(\mu_X, \sigma^2_X)$ and $Y \sim N(\mu_Y, \sigma^2_Y)$ ($X$ and $Y$ are independent), let $$Z = \begin{pmatrix} 2X + Y \\ X - 3Y \end{pmatrix}$$

I can easily calculate $\mu_{Z_{1,1}}, \mu_{Z_{2,1}}, \sigma^2_{Z_{1,1}}, \sigma^2_{Z_{2,1}}$, but I am having trouble figuring out how to calculate $\text{Cov}(Z_{1,1}, Z_{2,1})$ i.e. $\text{Cov}(2X + Y, X - 3Y)$. I understand that this is equal to $$E[(2X + Y)\times(X - 3Y)] - E[2X + Y]E[X - 3Y]$$ but am having trouble getting there, partially because I am worried that the two Gaussian random variables in $Z$ are functions of the original 2 Gaussian RVs $X, Y$.

Could somebody help me figure out how to find the covariance of this random vector?

$\operatorname{Cov}(X,Y)$ is linear in both $X$ and $Y$. So by linearity, $$\operatorname{Cov}(2X+Y,X-3Y) = 2\operatorname{Cov}(X,X)-5\operatorname{Cov}(X,Y) - 3\operatorname{Cov}(Y,Y) \\= 2\operatorname{Var}(X)-5\operatorname{Cov}(X,Y) - 3\operatorname{Var}(Y).$$ (By the way, computing the variances of $2X+Y$ and $X-3Y$ also requires using linearity.)