Joint density of $X_{(n)}-X_{(1)}$ and $X_{(n)}$ for Uniform distribution $[0,1]$

I am trying to find the Joint density of $X_{(n)}-X_{(1)}$ and $X_{(n)}$ for Uniform distribution $[0,1]$, where $X_{(1)}$ and $X_{(n)}$ are the smallest and largest observation (of $n$ samples).

Firstly, I found the density of $X_{(n)}, X_{(1)}$. Let $U=X_{(n)}$ and $V = X_{(1)}$.

\begin{aligned} f_{U,V}(u,v)=n(n-1)(u-v)^{n-2} \qquad ; 0<v<u<1 \end{aligned} Also the density of $X_{(n)}$ is

\begin{aligned} f_{U}(u)=nu^{n-1} \qquad ; 0<u<1 \end{aligned}

Now, I need to find the density of $R= X_{(n)}-X_{(1)}$.

Could you please give me an idea to proceed?

You can obtain the desired joint distribution by integrating like below. Note that for the sake of my sanity I have switched $U$ and $V$ :)
First, an aside for future readers who don't spot where your joint density $f_{U,V}$ comes from. Observe that $\mathbb{P}(U \ge u, V \le v) = \mathbb{P}(u \le X_i \le v$ for all $i) = (v-u)^n$ by independence of the variables. Now take the (negative) derivative in both variables.
Suppose $a < b$ (you will see how to get the full answer with a simple modification below). We have: $$\mathbb{P}(V - U \ge a, V \ge b) = \int_b^1 \mathbb{P}(U \le v - a, V \in dv)\\ = \int_b^1 \int_0^{v-a} n(n-1)(v-u)^{n-2} du dv \\ = -n\int_b^1 \left.(v-u)^{n-1}\right|_0^{v-a} dv \\ = (1-b^n) - n(1-b)a^{n-1},$$ where we used the condition $a < b$ to to ensure $v > b> a$ so that $v-a > 0$ in the integration limit in the second line. If $a > b$ you can split up the outer integral and see the part between $a$ and $b$ vanishes.
Towards the end of your post you refer only to $R = V-U$, the marginal. The density can be obtained by taking $b=0$ and then taking the (negative) derivative in $a$ after you include the case $a > b$ above. You find it is: $$n(n-1)(1-a)a^{n-2}$$
When $n=10$, this looks like the below. The histogram is from $50k$ simulations.
The next plot is the density for $n$ up to 10.