Find joint density function of X and X+Y (exponential distribution)
Solution 1:
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} $$
Solution 2:
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}$$