Even without computing anything, one can guess at the onset that $(Y,Z)$ is probably not independent since $Y\leqslant Z$ almost surely. Now, the most direct way to compute the distribution of $(Y,Z)$ might be to note that, for every $(y,z)$,
$$
[y\lt Y,Z\leqslant z]=\bigcap_{i=1}^n[y\lt X_i\leqslant z],
$$
hence, for every $0\lt y\lt z\lt1$,
$$
P(y\lt Y,Z\leqslant z)=P(y\lt X_1\leqslant z)^n=(z-y)^n.
$$
In particular, using this for $y=0$ yields, for every $z$ in $(0,1)$,
$$
F_Z(z)=P(Z\leqslant z)=z^n.
$$
Thus, for every $(y,z)$ in $(0,1)$,
$$
F_{Y,Z}(y,z)=P(Y\leqslant y,Z\leqslant z)
$$
is also
$$
F_{Y,Z}(y,z)=P(Z\leqslant z)-P(y\lt Y,Z\leqslant z)=z^n-(z-y)^n\cdot\mathbf 1_{y\lt z}.
$$
One may find more convenient to describe the distribution of $(Y,Z)$ by its PDF $f_{Y,Z}$, obtained as
$$
f_{Y,Z}=\frac{\partial^2F_{Y,Z}}{\partial y\partial z}.
$$
In the present case, for every $n\geqslant2$,
$$
f_{Y,Z}(y,z)=n(n-1)(z-y)^{n-2}\mathbf 1_{0\lt y\lt z\lt1}.
$$
For the first, $X-Y$ we need $x\in [0,\infty)$ and $(x-z)\in[0,1]$ where $z\in[-1,\infty)$
That's $-1\leq \max(0,z) \leq x\leq 1+z <\infty$, so the integration is: $$f_{X-Y}(z) = \mathbf 1_{z\in[-1,0)}\int\limits_0^{1+z}f_X(x)f_Y(x-z)\operatorname d x +\mathbf 1_{z\in[0,\infty)}\int\limits_z^{1+z}f_X(x)f_Y(x-z)\operatorname d x $$
For the second, $Y-X$, we need $x\in [0,\infty)$ and $x+z\in[0,1]$ where $z\in (-\infty, 1]$
That's $\max(0,-z) \leq x \leq 1-z < \infty$, so the integration is: $$f_{Y-X}(z) = \mathbf 1_{z\in[-\infty,0)}\int\limits_{-z}^{1-z}f_X(x)f_Y(x+z)\operatorname d x +\mathbf 1_{z\in[0,1]}\int\limits_0^{1-z}f_X(x)f_Y(x+z)\operatorname d x $$
Best Answer
Hint:
For fixed $z$:$$\mathsf P(Z<z)=\int_0^1\int_0^1[ax-y<z]dydx$$ where $[ax-y<z]$ is the function $\mathbb R^2\to\mathbb R$ taking value $1$ if $ax-y<z$ and taking value $0$ otherwise.