Find Distribution and Conditional Expectation / Variance of Multivariate Gaussian random variables I have that
$$
X \in N \begin{pmatrix} \begin{pmatrix} 1 \\ 1 \end{pmatrix},\begin{pmatrix} 3 & 1 \\ 1 & 2 \end{pmatrix} \end{pmatrix}
$$
and would like to compute the distribution of $X_1 + X_2 | X_1-X_2 $.
I have started by setting,
$$
Y_1=X_1+X_2 \\Y_2=X_1-X_2
$$ 
and thus get that,
$$
\mu_1=E[Y_1]=E[X_1]-E[X_2]=2 \\ \mu_2=E[Y_2]=E[X_1]-E[X_2]=0.
$$
I know that for the case $n=2$ the Conditional Expectation and Variance can be computed using,
$$ 
E[Y_1|Y_2]=\mu_1 + \frac{\sigma_{12}}{\sigma_2^2}(y_2-\mu_2)\\
Var(Y_1|Y_2) =\sigma_1^2-\frac{\sigma_{12}^2}{\sigma_2^2}
$$
So to compute variances $\sigma_1^2$ and $\sigma_2^2$, I start by setting,
$$ 
\sigma_1^2=Var(X_1+X_2)=E[(X_1+X_2)^2]-(E[X_1-X_2])^2=E[X_1^2]+2E[X_1]E[X_2]+E[X_2^2]-(E[X_1]-E[X_2])^2
$$
But to compute the second moment of $X$, i.e. $E[X^2]$, I want to use that,
$$
E[X^n]=\frac{1}{i^2}\frac{d^n}{dt^n}\varphi_X(t=0),
$$
where the characteristic function for a Normal random variable is $\varphi_X(t)=e^{i\mu t-\frac{1}{2}t^2 \sigma^2}$. However I get that, 
$$
E[X^2]=\frac{1}{i^2}\frac{d^2}{dt^2}e^{i\mu 0-\frac{1}{2}0 \sigma^2}=0.
$$
This would then imply that,
$$ 
\sigma_1^2=Var(X_1+X_2)=0+2\cdot1\cdot1+0-(1-1)^2=2\neq 7.
$$
This does not correspond to the answer which states that,
$$ 
E[Y_1|Y_2]=2 + \frac{1}{3}(y_2-0)\\
Var(Y_1|Y_2) =7-\frac{1}{3}=\frac{20}{3}
$$
Could anybody please show me the correct way and explain the meaning of writing (i.e. what does B do and how can I use it?) as done here,
$$
Y=\begin{pmatrix} Y_1 \\ Y_2  \end{pmatrix}=\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}X=BX.
$$
Tanks!
 A: You have a few mistakes in your analysis.
First, your alarm should sound when you obtain a result like $E[X^2]=0$ ($X$ is a vector; did you mean $X_1$?).
Having $E[X_1^2]=0$ would imply that $X_1=0$ almost surely, which it's not.
As a side note, you can compute $E[X_1^2]$ fairly easily using
$$E[X_1^2]=Var(X_1)+E[X_1]^2,$$
without resorting to the characteristic function.
Second, in your calculation of the variance of $X_1+X_2$, you assume that
$$
E[X_1X_2]=E[X_1]E[X_2].
$$
This is not true in general, rather only when the two variables are uncorrelated. For instance, if $X_1=+1$ or $-1$ with probability $1/2$, and $X_2=-X_1$, then $E[X_1]E[X_2]=0$ but $E[X_1X_2]=-1$.

What you're dealing with here is an affine transformation of a multivariate normal variable.
Matrices are your friend; use them.
Writing $Y=BX$ makes all your analysis much simpler and straightforward.
If $X\sim\mathcal{N}(\mu,\Sigma)$, then the distribution of $Y=c+BX$ is:
$$
Y \sim \mathcal{N}\left( c+B\mu, B\Sigma B^T \right).
$$
In your case, you can stack $Y_1$ and $Y_2$ into one vector $Y=(Y_1,Y_2)^T$, giving
$$
Y = \begin{pmatrix}1&1\\1&-1\end{pmatrix}X = BX \sim\mathcal{N}\left(
B\begin{pmatrix}1\\1\end{pmatrix},
B\begin{pmatrix}3&1\\1&2\end{pmatrix}B^T
\right).
$$
Thus the covariance matrix of $Y$ is easily calculated to be
$$
B\begin{pmatrix}3&1\\1&2\end{pmatrix}B^T
= \begin{pmatrix}1&1\\1&-1\end{pmatrix}\begin{pmatrix}3&1\\1&2\end{pmatrix}\begin{pmatrix}1&1\\1&-1\end{pmatrix}
= \begin{pmatrix}1&1\\1&-1\end{pmatrix}\begin{pmatrix}4&2\\3&-1\end{pmatrix}
= \begin{pmatrix}7&1\\1&3\end{pmatrix}.
$$
You should then be able to get the correct result.
