Assuming that $X$ and $Y$ are independent, the joint distribution is
$f(x,y) = \frac{1}{\theta_1 \theta_2}$ if $0 \le x \le \theta_1, 0 \le y \le \theta_2$ and $0$ otherwise.
The correct way to do this is to integrate $f(x,y)$ over the region $x-y\le z$. The easiest way to compute this integral is to do it geometrically. Compute the area of the intersection of the region $x-y\le z$ and the rectangle $[0, \theta_1] \times [0, \theta_2]$ and divide it by $\theta_1 \theta_2$.
You can also use the convolution to compute the distribution of $Z = X - Y$ and find $P\{Z \le z\}$ by integrating $\int_0^z f_Z(s) ds$.
You will need to draw at least one picture. First we identify the part $L$ of the plane where the joint density function of $X$ and $Y$ "lives."
Draw the vertical lines $x=1$ and $x=2$. Draw the line $y=x$. The joint density lives on the part $L$ of the first quadrant that is between $x=1$ and $x=2$, and below $y=x$. This region $L$ is a trapezoid with corners $(1,0)$, $(2,0)$, $(2,2)$, and $(1,1)$.
Let $Z=Y-X$. We want the cdf $F_Z(z)$ of $Z$, that is, $\Pr(Z\le z)$. So we want $\Pr(Y\le X+z)$. We want the probability that $(X,Y)$ lands in the part $D$ of $L$ that lies below the line $y=x+w$.
For fixed $z$, we will draw the line $y=x+z$. The only interesting values of $z$ are between $-2$ and $0$. Note that the line $y=x+z$ is always parallel to the line $y=x$.
The geometry is different for small negative $z$ than for large negative $z$. For if $z$ is small negative, the region $D$ below $y=x+z$ but in $L$ is a trapezoid.
At the value of $z$ for which the line $y=x+z$ passes through $(1,0)$, the part $D$ of $L$ below $y=x+z$ turns into a triangle. This happens when $z=-1$.
Now we can calculate. For $-1\lt z\le 0$, we integrate first with respect to $y$, $y=0$ to $y=x+z$, and then with respect to $x$, from $1$ to $2$. So we have
$$F_Z(z)=\int_{x=1}^2\left(\int_{y=0}^{x+z}\frac{3}{7}x\,dy \right)\,dx$$
when $-1\lt z\le 0$.
For $-2\le z\le -1$, we are integrating over a triangle. Note that the triangle has corners $(-z,0)$, $(2,0)$, and $(2,2+z)$.
So as before, in the inner integration, $y$ will go from $0$ to $x+z$. Then for the outer integration, $x$ goes from $-z$ to $2$. So
$$F_Z(z)=\int_{x=-z}^2\left(\int_{y=0}^{x+z}\frac{3}{7}x\,dy \right)\,dx$$
when $-2\le z\le -1$.
The rest is mechanical. Compute the integrals, and differentiate to find $f_Z(z)$.
Remark: There are other ways to solve the problem. One can find the density function of $Z$ directly, by using a variant of the usual convolution procedure. Doing this still requires careful attention to the geometry, to get the limits right.
Best Answer
We have $$\begin{align*} F_Z(z) &= P(X+Y \leq z) = E[\mathbf{1}_{X+Y \leq z}]\\ &= \int_{[0,\infty) \times [0,\infty)} \mathbf{1}_{x+y \leq z} \; f_{X,Y}(x,y) \, \mathrm{d}x \mathrm{d}y \end{align*}$$ using the following property of the density: $E[g(X,Y)] = \int_{[0,\infty) \times [0,\infty)} g(x,y)\, f_{X,Y}(x,y) \, \mathrm{d}x \mathrm{d}y$. Now notice that $\mathbf{1}_{x+y \leq z} = \mathbf{1}_{x \leq z} \mathbf{1}_{y \leq z - x}$, that is $x+y \leq z$ if and only if $x \leq z$ and $y \leq z - x$. So using Fubini's theorem, we get $$ \begin{align*} F_Z(z) = \int_{[0,\infty) \times [0,\infty)} \mathbf{1}_{x \leq z} \mathbf{1}_{y \leq z - x} \; f_{X,Y}(x,y) \, \mathrm{d}x \mathrm{d}y &= \int_{[0,\infty)} \mathbf{1}_{x \leq z} \left(\int_{[0,\infty)}\mathbf{1}_{y \leq z - x} \; f_{X,Y}(x,y) \, \mathrm{d}y \right)\mathrm{d}x \\ &= \int_0^z \left(\int_0^{z-x} f_{X,Y}(x,y) \, \mathrm{d}y \right)\mathrm{d}x \end{align*}$$ which proves the result. Normally you would not have to go into as much detail and from the first equation above you could directly obtain the result.