Let $X\sim N(0,1)$ and $Y\sim N(1,4)$ two random independent normal variables.
We have to find the joint density of $\begin{pmatrix}X-Y \cr Y\end{pmatrix}$.
I defined a diffeomorphism $\varphi$ such that
$$\begin{pmatrix}X-Y \cr Y\end{pmatrix}=\varphi\begin{pmatrix}X \cr Y\end{pmatrix}$$
as
$$\varphi:\mathbb{R}\times \mathbb{R}\to \mathbb{R}\times \mathbb{R}, \quad \varphi(x,y)=(x-y,y)$$
and I determinated the inverse $\varphi^{-1}$: given $(s,t)\in \mathbb{R}\times\mathbb{R}$, I searched $(x,y)\in \mathbb{R}\times\mathbb{R}$ such that
$$\begin{cases}s=x-y \cr t=y\end{cases} \Longrightarrow \begin{cases}x=s+t \cr y=t\end{cases}$$
So, $\varphi^{-1}(s,t)=(s+t,t)$ and its jacobian matrix is
$$ J_{\varphi^{-1}}(s,t)=\begin{pmatrix}1 & 0 \cr 1 & 1\end{pmatrix}$$
which has determinant equals to 1.
So, If I say $S=X-Y$ and $T=Y$, I have that
$$f_{(S,T)}(s,t)=f_{(X,Y)}(\varphi^{-1}(s,t))\cdot |\mathrm{det} J_{\varphi^{-1}}(s,t)|=f_{(X,Y)}(s+t,t)$$.
Now I used the fact that $X$ and $Y$ are independent and I can write that
$$f_{(X,Y)}(s+t,t)=f_X(s+t)f_Y(t)$$
I have that
$$f_X(s+t)=\frac{1}{\sqrt{2\pi}}e^{-\frac{(s+t)^2}{2}}$$
and
$$f_Y(t)=\frac{1}{\sqrt{8\pi}}e^{-\frac{(t-1)^2}{8}}$$
Then, I obtained that
$$f_{(S,T)}(s,t)=\frac{1}{\sqrt{2\pi}}e^{-\frac{(s+t)^2}{2}}\frac{1}{\sqrt{8\pi}}e^{-\frac{(t-1)^2}{8}}$$
but I can't rewrite this as a notable density, suggestions??
It is known that if $X,Y$ are jointly normal, then any linear transformation is also jointly normal. As suggested in the comments and answers, it is better to do this in a matrix way. You can use the following: If you are interested in a transformation $A\textbf{X}$, then this transformation must be normal with a mean vector $\boldsymbol{\mu}^{*} = A \boldsymbol{\mu}$ and a covariance matrix $\boldsymbol{\Sigma}^{*} = A Var(\textbf{X}) A^T$. In our case $$ \begin{pmatrix} X_1 \\\ X_2 \end{pmatrix} \sim \mathcal{N}_2\left( \begin{pmatrix} 0 \\\ 1 \end{pmatrix},\begin{pmatrix} 1 & 0 \\\ 0 & 4 \end{pmatrix} \right) $$ we are interested in the transformation $$ \begin{pmatrix} X_1 & X_2 \\\ 0 & -X_2 \end{pmatrix} = \begin{pmatrix} 1 & -1 \\\ 0 & 1 \end{pmatrix} \begin{pmatrix} X_1 \\\ X_2 \end{pmatrix} = A \textbf{X} $$
now we only calculate what is necessary $$ \boldsymbol{\mu}^{*} = A \boldsymbol{\mu} = \begin{pmatrix} 1 & -1 \\\ 0 & 1 \end{pmatrix} \begin{pmatrix} 0 \\\ 1 \end{pmatrix} = \begin{pmatrix} -1 \\\ 1 \end{pmatrix} $$
$$ \boldsymbol{\Sigma}^{*} = A Var(\textbf{X}) A^T =\begin{pmatrix} 1 & -1 \\\ 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\\ 0 & 4 \end{pmatrix} \begin{pmatrix} 1 & 0 \\\ -1 & 1 \end{pmatrix} = \left(\begin{matrix} 5 & -4 \\ -4 & 4 \end{matrix}\right) $$ finally we have
$$ \begin{pmatrix} X_1 - X_2 \\\ X_2 \end{pmatrix} \sim \mathcal{N}_2 \left(\begin{pmatrix} -1 \\\ 1 \end{pmatrix} , \left(\begin{matrix} 5 & -4 \\ -4 & 4 \end{matrix}\right)\right) $$ or in your notation $$ \begin{pmatrix} X- Y \\\ Y \end{pmatrix} \sim \mathcal{N}_2 \left(\begin{pmatrix} -1 \\\ 1 \end{pmatrix} , \left(\begin{matrix} 5 & -4 \\ -4 & 4 \end{matrix}\right)\right) $$ Now, you were wondering how to write the density found in a familiar way, well, use the Wikipedia (Henry`s comment) page to put together the form of the bivariate normal with the parameters obtained. It is complex, but not impossible. If you are interested in recognizing the distribution in this type of situation, it is better to proceed in a matrix manner.