Find joint density function of X and X+Y (exponential distribution) I tried finding the following density function:

Let $X$ and $Y$ be independent random variables, each having the exponential distribution with parameter $\lambda$. Find the joint density function of $X$ and $Z=X+Y$.

$\begin{aligned}f_{X,X+Y}=&\frac{\partial}{\partial x\partial z}\mathbb P (X\leq x,Y\leq z-x)\\=&\int_{-\infty}^x\int_{-\infty}^{z-x}f(u,v)\,\mathrm dv=\frac{\partial}{\partial z}\left[\frac{\partial}{\partial x}\int_0^x\int_0^{z-x}\lambda^2e^{-\lambda u-\lambda v}\,\mathrm dv\,\mathrm d u\right]\\
=&\frac{\partial}{\partial z}\int_0^{z-x}\lambda^2e^{-\lambda x-\lambda v}\,\mathrm dv=\frac{\partial}{\partial z}\int_x^{z}\lambda^2e^{-\lambda x-\lambda(v-x)}\,\mathrm dv\\
=&\lambda^2e^{-\lambda x-\lambda(z-x)}=\lambda^2e^{-\lambda z}
\end{aligned}$
But I don't think this is correct, because $\mathbb P(X\leq x,X+Y\leq z)$ isn't equal to $\mathbb P(X\leq x,Y\leq z-x)$. But I can't think of anything else.
Could someone help me a little bit with this exercise?
I'm going to try the following:
$\mathbb P(X\leq x,Z\leq z)=\int_{-\infty}^x\int_{-\infty}^zf_X(u)f_Z(v)\,\mathrm du\,\mathrm dv$,
where I use the convolution formula. If I made no mistake, this gives
$f_{X+Y}(z)=z\lambda^2e^{-\lambda z}$, if $z>0$.
Then calculating the joint density function:
$\begin{aligned}f_{X,X+Y}(x,z)=\frac{\partial}{\partial z}\left[\int_{-\infty}^z\lambda^3ve^{-\lambda x-\lambda v}\,\mathrm dv\right].
\end{aligned}$
Oh wait, this is also incorrect, because $X$ and $X+Y$ are dependent.
I really don't know how to proceed; could someone help?
I've also looked at this post, but they mention the Jacobian, and my teacher has chosen to skip the paragraph about the Jacobian, so it wouldn't make sense if I had to use that.
 A: Directly, by the Jacobian change of variables, we have:
$$\begin{align}f_{X,X+Y}(x,z) ~&=~ f_{X,Y}(x,z-x)~\lVert\tfrac{\partial(x,z-x)}{\partial (x,z)}\rVert\\[1ex]& =~ \lambda^2\mathsf e^{-\lambda(x+z-x)}\mathbf 1_{(x,z)\in\Bbb R^2\wedge0\leqslant x\leqslant z}\\[1ex] &=~ \lambda^2\mathsf e^{-\lambda z}\mathbf 1_{(x,z)\in\Bbb R^2\wedge 0\leqslant x\leqslant z}\end{align}$$
(The Jacobian determinant is, conveniently, $1$ in this case).

By your method. (Using the Fundamental Theorem of Calculus:)
$$\begin{align}f_{X,X+Y}(x,z) ~&=~ \dfrac{\partial^2\qquad }{\partial~x~\partial~z}\mathsf P(X\leq x,X+Y\leq z) \\[1ex] &=~ \dfrac{\partial^2\qquad }{\partial~x~\partial~z} \int_0^x f_X(s)\mathsf P(Y\leq z-X\mid X=s)\operatorname d s \\[1ex] &=~ f_X(x)\,\dfrac{\partial\quad}{\partial~ z~}\mathsf P(Y\leq z-x) \\[1ex] &= f_X(x)\,\dfrac{\partial\quad}{\partial~ z~}\int_0^{z-x}f_Y(t)\operatorname d t \\[1ex] &=~ f_X(x)\,f_Y(z-x) \\[1ex] &=~ \text{(see above)}\end{align} $$
A: In this solution I try to do everything step by step without referring to more advanced knowledge. The most honorable OP has to accept that this problem is not that simple. (Elementary but not simple; not tricky but needs work.)
Let's concentrate on the common distribution function:
$$F_{X,X+Y}(x,y)=P(X<x\cap X+Y<y)=$$
$$=\int_0^{\infty}P(X<x\ \cap\ X+Y<y \mid X=z)\ \lambda e^{-\lambda z} \ dz.\tag1$$
After substituting $z$ for $X$ we get
$$P(X<x\ \cap\ X+Y<y \mid X=z)=P(z<x\cap Y<y-z)=\mathbb I_{z<x}P(Y<y-z)$$
because  $X$ and $Y$ are independent.
Since $\mathbb I_{z<x}P(Y<y-z)=0$ if $\min(x,y)\leq z$, $(1)$ will be
$$\lambda\int_0^{\min(x,y)}P(Y<y-z)  e^{-\lambda z}\ dz=$$
$$=\lambda\int_0^{\min(x,y)}\left(1-e^{-\lambda(y-z)}\right)e^{-\lambda z}\ dz=$$
$$=\lambda\int_0^{\min(x,y)}e^{-\lambda z}\ dz-\lambda e^{-\lambda y}\int_0^{\min(x,y)}\ dz=$$
$$=1-e^{-\lambda \min(x,y)}-\lambda e^{-\lambda y}\min(x,y).$$
Now, let's compute the common density. For $x<y$ we have 
$$f_{X,X+Y}(x,y)=\frac{\partial^2 F_{X,X+Y}}{\partial y \partial x}=\frac{\partial}{\partial y}\frac{\partial}{\partial x}\left[1-e^{-\lambda x}-\lambda e^{-\lambda y}x\right]=$$
$$=\lambda\frac{\partial}{\partial y}\left(e^{-\lambda x}-e^{-\lambda y}\right)=\lambda^2e^{-\lambda y}.$$
But if $y\leq x$ then the operation above results in $0$.
As a result we have 
$$f_{X,X+Y}(x,y)=\begin{cases}
\lambda^2e^{-\lambda y}&\text{ if } 0\leq x\leq y\\
0& \text{ otherwise}.
\end{cases}$$
