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.)
