This problem statement has been featured in a different question, but the user who posted that question was asking about clarification of the definitions within the problem statement, while I am trying to verify my solution, so I don't believe this should be considered as a repost. One important nuance here is that we're supposed to solve this problem without having been 'taught' anything about order statistics and sample sizes, just treating it as a problem like any other joint distribution problems that we have been doing (i.e. just figuring out the region of integration, and solving the integrals.)
My attempt:
Define $R = |X_2 - X_1|$.
For $r\in(0,1)$, $P(R<r) = P(|X_2 - X_1| < r) = 2P(X_2 - X_1 < r)$ by symmetry.
Computing $P(X_2 - X_1 < r):$
$\displaystyle{P(X_2 - X_1 < r) = \iint_{x_2-x_1<r}f_{X_1}(x_1)f_{X_2}(x_2)dx_1dx_2 = \int_0^{r}\int_{0}^{r+x_1}(2x_1)(2x_2)dx_2dx_1}$.
I am a little unconfident that I chose the correct limits of integration, and that I set up everything else correctly so far. But if I did, then computing this integral yields $\frac{17}6r^4, r\in(0,1)$. Thus, (and apologies for my poor MathJax):
$F_R(r) = 1$ when $r\geq1$,
$F_R(r) = \dfrac{17}3r^4$ when $r\in(0,1)$,
$F_R(r) = 0$ when $r\leq0$.
Differentiation yields that:
$f_R(r) = \dfrac{68}3r^3$ when $r\in(0,1)$,
$f_R(r) = 0$ when $r\notin(0,1)$.
Everything check out?
We are interested in the probability mass in the region $R = \{(x_1,x_2): |x_1 - x_2| \leq r, x_1 \in [0,1], |x_2| \in [0,1] \}$. This is essentially the region in the unit square between the lines $x_2 = x_1 - r$ and $x_2 = x_1 + r$. To simplify our calculation, we will calculate the region in the unit square outside $R$. Since, by symmetry the mass of the region below the line $x_2 = x_1 - r$ is the same as that of above the line $x_2 = x_1 + r$ (as you have already noted), we will only compute the former. Denoting it by $p$, we have, \begin{align*} p & = \int_{r}^1 \int_{0}^{x_1 - r} (2x_1) \cdot (2x_2) \ \mathrm{d}x_2 \ \mathrm{d}x_1 \\ & = \int_{r}^1 2x_1(x_1 - r)^2 \ \mathrm{d}x_1 \\ & = 2\int_{r}^1 [r(x_1 - r)^2 + (x_1 - r)^3] \ \mathrm{d}x_1 \\ & = \frac{2r(1 - r)^3}{3} + \frac{(1 - r)^4}{2}. \end{align*}
Thus, the required probability is $1 - 2p$.