Assuming $$ X\sim N_2(\mu_1, \sigma_1)$$ and $$Y \sim N_2(\mu_2, \sigma_2) $$
Would the joint distribution of
\begin{bmatrix} X -Y \\ X+Y \end{bmatrix}
be 4 dimensional? It would be multivariate normal too because it is the joint of 2 Multivariate normal vectors, correct?
No. Unless your notation is misleading, $X$ and $Y$ are univariate.
Consequently each of $D=X-Y$ and $S=X+Y$ are univariate.
A vector made of those two random variables $[D,S]^\top$ is bivariate.
There are too many incorrect notions in your assertions. Regardless of whether $X$ and $Y$ are multivariate normal or univariate normal random variables, $X+Y$ and $X-Y$ are not necessarily normal random variables (multivariate or univariate as the case may be). If
and
then $\left[\begin{array}{c}X+Y\\X-Y\end{array}\right]$ is a $2n$-variate normal random variable. Notice that $\left[\begin{array}{c}X+Y\\X-Y\end{array}\right]$ is a linear transformation of $\left[\begin{array}{c}X\\Y\end{array}\right]$, and it is a defining property of multivariate normality that linear transformations of multivariate normal random variables result in multivariate normal random variables. To forestall the inevitable question "Didn't I say that $X$ and $Y$ are multivariate normal and doesn't that automatically mean that $\left[\begin{array}{c}X\\Y\end{array}\right]$ is also multivariate normal?", note that $$\left[\begin{array}{c}X\\Y\end{array}\right] 2n\text{-variate normal} \implies X, Y ~ n\text{-variate normal}$$ but the reverse implication is not true: $$\require{cancel}X, Y ~ n\text{-variate normal} \cancel{\implies} \left[\begin{array}{c}X\\Y\end{array}\right] 2n\text{-variate normal}$$ Many counterexamples for the alleged reverse implication can be found in this answer by Moderator cardinal to the question Is it possible to have a pair of Gaussian random variables for which the joint distribution is not Gaussian?
