Find the pdf of a bivariate transformation using change of variables Problem
$X,Y\sim n(0,1)$. $X,Y$ are indpendent.  $U = X+Y$, $V=X-Y$. Find $f_{U,V}(u,v).$ 

Work
I've made some progress, but there are a number of steps remaining and I'm not sure how to proceed. I know there's a formula for this sort of thing 
$$f_{U,V}(u,v) = f_{X,Y}(h1(u,v), h2(u,v))|J|$$
but I wanted to try working it out as explicitly as possible, using the definitions and integral change of variables in order to see what's going on under the hood.
$$f_{U,V}(u,v)$$
$$= \frac{\partial}{\partial u \partial v} F_{U,V}(u,v)$$
by the second FTC.
$$= \frac{\partial}{\partial u \partial v} \int_{\{(u,v): U < u, V < v\}} f_{U,V}(u,v)dudv $$ by the definition of $F_{U,V}$.
So we'll change variables from U-V into X-Y coordinates. First we'll measure the area element
$$dudv = \frac{\partial(u,v)}{\partial(x,y)}dxdy=2dxdy$$ using the Jacobian determinant, so that we can rescale appropriately the sum of the various volume chunks.
Second, we change the integrand. Do we plug in values of $u=x+y$ and $v=x-y$ in the integrand? This seems like it would make sense given that we want to replace $u$ and $v$ with $x$ and $y$, and this is the correct equivalence between the two pairs.
Third, we change the region. 
$$R=\{(u,v):U < u, V < v\}$$
$$=\{(u,v):X+Y < u, X-Y < v\}$$
but I'm not sure what we're getting out of this.
Fourth, we have to do something with the partial derivatives. Perhaps chain rule, or something cool with the FTC -- but I don't know how that would work if we've already changed out our $u$ and $v$ that we're supposed to be differentiating with respect to.

A big part of learning how to do this is surely in getting used to the notation and the manipulations, so I'd most appreciate as explicit a solution as possible, showing how each manipulation works and is motivated.
 A: We know $$\begin{align} F_{U,V}(u,v)&=\mathbb{P}(U<u,V<v) \\ &=\mathbb{P}(X+Y<u,X-Y<v)\\ &=\frac {1}{2\pi}\iint_Re^{-\frac 12 (x^2+y^2)}\mathrm{d}x\mathrm{d}y\end{align}$$
Where $$\begin{align} R&=\{(x,y):-\infty<x,y<\infty, x+y<u, x-y<v\}\\ &= \{(x,y):-\infty<y<\infty, -\infty<x<\min(u-y,v+y)\}  \\ &=\{(x,y):-\infty<y<\infty, -\infty<x<u-y,u-y<v+y\}\cup \{(x,y):-\infty<y<\infty, -\infty<x<v+y,v+y<u-y\} \\ &=\{(x,y):-\infty<x<u-y,\frac{u-v}{2}<y<\infty\}\cup \{(x,y): -\infty<x<v+y,-\infty<y<\frac{u-v}{2}\} \end{align}$$
We thus get an expression for the integral as
$$F_{U,V}(u,v)=\frac{1}{2\pi}\left(\int_{\frac{u-v}{2}}^{\infty}\int_{-\infty}^{u-y}e^{-\frac 12 (x^2+y^2)}\mathrm{d}x\mathrm{d}y+\int_{-\infty}^{\frac{u-v}{2}}\int_{-\infty}^{v+y}e^{-\frac 12 (x^2+y^2)}\mathrm{d}x\mathrm{d}y\right)$$
Then take partial derivatives to obtain $f_{U,V}$. As an example of how to obtain the partial with respect to $u$ of the first term, set $$F(u,y)=\int_{-\infty}^{u-y}e^{-\frac 12 (x^2+y^2)}\mathrm{d}x$$ and hence $$\frac{\partial}{\partial u}\int_{\frac{u-v}{2}}^{\infty}F(u,y)\mathrm{d}y=\int_{\frac{u-v}{2}}^{\infty}\frac{\partial}{\partial u}F(u,y)\mathrm{d}y-\frac 12 F(u,\frac{u-v}{2})$$
and 
$$\frac{\partial}{\partial u}F(u,y)=\frac{\partial}{\partial u}\int_{-\infty}^{u-y}e^{-\frac 12 (x^2+y^2)}\mathrm{d}x=e^{-\frac 12 ((u-y)^2+y^2)}$$
By the Leibniz rule and FTC respectively. Now do this for the second term, and then repeat this process with respect to $v$. 
A: Comment (continued). Here is a simulation with 50,000 points, illustrating the rotation, reflection,
and stretching when $(X,Y)$ is transformed to $(U, V).$ Each point
has the same color before and after transformation.
par(mfrow=c(1,2), pty="s")
  x = rnorm(50000);  y = rnorm(50000)
  u = x + y;  v = x - y
  b = max(-min(x,y), max(x,y))
  cx = (x + b); cy = (y+b)  # colors
  plot(x,y, xlim=c(-6,6), ylim=c(-6,6), col=rgb(cx, b, cy, max=2*b), pch=".")
    abline(v = 0, col="darkgreen"); abline(h = 0, col="darkgreen") # axes
  plot(u,v, xlim=c(-6,6), ylim=c(-6,6), col=rgb(cx, b, cy, max=2*b), pch=".")
    abline(a = 0, b=-1, col="darkgreen") # rotated
    abline(a = 0, b=1, col="darkgreen")  #   axes
par(mfrow=c(1,1), pty="m")


Outline of PDF transformation method (added later): The Jacobian of the inverse transformation is $J = -1/2,$ so
$$f_{(U,V)}(u,v) = \frac{|J|}{2\pi}\exp\left[-\frac{1}{2}\left(\frac{u+v}{2}\right)^2  -\frac{1}{2}\left(\frac{u-v}{2}\right)^2\right] = \cdots\\
= \frac{e^{-u^2/4}}{\sqrt{2}\sqrt{2\pi}} \times \frac{e^{-v^2/4}}{\sqrt{2}\sqrt{2\pi}}.$$
This implies that $U$ and $V$ are independent and each is distributed $Norm(\mu=0, \sigma=\sqrt{2}).$ (Acknowledgment: This is mostly copied from
Wackerly et al: Mathematical Statistics with Applications, 
7e, Cengage, 2008, p326.)  
