If $X$ and $Y$ have joint density function $$ f(x,y)=\begin{cases} \frac{1}{2}(x+y)e^{-x-y}&\text{if }x,y>0,\\ 0&\text{otherwise,} \end{cases} $$ find the density function of $Z=X+Y$.
I did the following: $$ \begin{aligned} f_Z(z)&=1/2\frac{\partial}{\partial z}\int_0^z\int_0^{z-x}(x+y)e^{-x-y}dydx\\ &=1/2\frac{\partial}{\partial z}\int_0^z\left[[-(x+y)e^{-x-y}]_0^{z-x}+\int_o^{z-x}e^{-x-y}\right]dx\\ &=1/2\frac{\partial}{\partial z}\int_0^z\left[-ze^{-z}+xe^{-x}+[-e^{-x-y}]_0^{z-x}\right]dx\\ &=1/2\frac{\partial}{\partial z}\int_0^z-ze^{-z}+xe^{-x}-e^{-z}+e^{-x}dx\\ &=1/2\frac{\partial}{\partial z}\left[\int_0^z-ze^{-z}-e^{-z}+e^{-x})+\int_0^z xe^{-x}dx \right]\\ &=1/2\frac{\partial}{\partial z}\left[(-z^2e^{-z}-ze^{-z}-e^{-z}+1)+\int_0^z xe^{-x}dx \right]\\ &=1/2\frac{\partial}{\partial z}\left[-z^2e^{-z}-ze^{-z}-e^{-z}+1-ze^{-z}-e^{-z}+1\right]\\ &=1/2\frac{\partial}{\partial z}\left[-z^2e^{-z}-2ze^{-z}-2e^{-z}+2\right]\\ &=1/2[-2ze^{-z}+z^2e^{-z}-2e^{-z}+2ze^{-z}+2e^{-z}]\\ &=1/2z^2e^{-z}. \end{aligned} $$
EDIT
Apparently, I can't use the fundamental theorem, but if I just finish the outer integral and then differentiate, then things go well.
Is there an easier approach though?
Here I present you another way of doing it, which is much faster.
Probability fact:
For $Z=X+Y$, and $D$ is the positive $x$-$y$ plane, \begin{align} f_Z(z) &=\int\int_D f(x,y)\delta(z-x-y) dxdy \\ & = \frac{1}{2}z e^{-z} \int\frac{1}{\sqrt{2}} dl \\ & =\frac{1}{2}z^2 e^{-z} \end{align}
In the second equality, we have used $\int _{\mathbf {R} ^{2}}f(\mathbf {x} )\,\delta (g(\mathbf {x} ))\,d\mathbf {x} =\int _{g^{-1}(0)}{\frac {f(\mathbf {x} )}{|\mathbf {\nabla } g|}}\,dl$, here $g(\mathbf{x})=z-x-y,$ and $dl$ is the line element where $g(\mathbf{x})=0$. The facor of $1/\sqrt{2}$ comes from $|\mathbf {\nabla } g|$. In the last equality, we simply know that the integral is the length of the line connecting $(z,0)$, $(0,z)$.