Using the following change of coordinates,
u(x, y, z) = y - x
v(x, y, z) = z - y
w(x, y, z) = z
we solve for x, y, and z in terms of u, v, and w :
x(u, v, w) = - u - v + w
y(u, v, w) = - v + w
z(u, v, w) = w
We then compute the Jacobian for the transformation,
[tex]J = \begin{bmatrix} x_u & x_v & x_w \\ y_u & y_v & y_w \\ z_u & z_v & z_w\end{bmatrix} = \begin{bmatrix} - 1 & -1 & 1 \\ 0 & -1 & 1 \\ 0 & 0 & 1 \end{bmatrix}[/tex]
so that det(J) = 1.
Now, the region D transforms to D', where
[tex]D' = \left\{ (u, v, w) : 0 \le u \le 2 \text{ and } 0 \le v \le 1 \text{ and } 0 \le w \le 4\right\}[/tex]
so the integral is
[tex]\displaystyle \iiint_D xy \, dV = \iiint_{D'} (-u-v+w) (-v+w) \det(J) \, dV \\\\ = \int_0^2 \int_0^1 \int_0^4 (uv - uw + v^2 - 2vw + w^2) \, dw \, dv \, du \\\\ = \int_0^2 \int_0^1 \left(4uv - 8u + 4v^2 - 16v + \frac{64}3\right) \, dv \, du \\\\ = \int_0^2 \left(\frac{44}3 - 6u\right) \, du \\\\ = \boxed{\frac{52}3}[/tex]