Analytically computing a double integral? I have used Mathematica to show that
$$\int_{0}^\infty\int_{0}^y x^{c-1}(y-x)^{L-c}e^{bx-(a+b+d)y}~dx~dy = \dfrac{(c-1)!(L-c)!}{(a + d)^c (a + b + d)^{L+1-c}}$$
where $L \ge c$ with $c,L$ positive integers and $a,b,d$ positive reals.
I need to show this by hand. Can someone point me in the right direction? References are welcome.
Thanks in advance for any help.
 A: Make the change of variable $\cases{u=x\\v=y-x} \iff  \cases{x=u\\y=u+v}$ (with jacobian 1).
In this way
$$\int_{y=0}^{\infty}\int_{x=0}^y x^{c-1}(y-x)^{L-c}e^{bx-(a+b+d)y}~dx~dy$$
becomes
$$\int_{v=0}^{\infty}\int_{u=0}^{\infty} u^{c-1}v^{L-c}e^{-bv}e^{-(a+d)(u+v)}~du~dv$$
giving the separated form:
$$\int_{u=0}^{\infty} u^{c-1}e^{-(a+d)u}~du\int_{v=0}^{\infty}v^{L-c}e^{-(a+b+d)v}~dv$$
each of these integrals being obtainable as a value of $\Gamma$ function.
Edit: justification of the domain $[0,+\infty)^2$ for  $(u,v)$. 
See figure below.
Constraint $0 \leq x \leq y \ \iff \ (x,y) \in T$ where $T$ is the non bounded triangular domain depicted in blue. 
Integration limits in $\int_{y=0}^{\infty}\int_{x=0}^y...$ mean that, for a given $y$, point $(x,y)$, "sweeps" an horizontal line segment (figure on the left). Turning to new variables $(u,v)$ generates, one has the sweeping mode represented on the right, along half lines with equations $y=u+v$, i.e., $y=x+v$ ; this explains that $y$ and $v$ (the so-called $y$ intercept) should begin at $0$. Said otherwise, giving any ($\geq 0$) value of $u$ (i.e., of $x$) and any ($\geq 0$) value of  $v$, we have a unique point in the triangular domain.

