Let $X_{(1)},X_{(2)},\ldots,X_{(n)}$ be the order statistics corresponding to $X_1,X_2,\ldots,X_n$.
Making the transformation $(X_{(1)},\ldots,X_{(n)})\mapsto (Y_1,\ldots,Y_n)$ where $Y_1=X_{(1)}$ and $Y_i=X_{(i)}-X_{(i-1)}$ for $i=2,3,\ldots,n$, we have $Y_i$ exponential with mean $1/(n-i+1)$ independently for all $i=1,\ldots,n$.
Therefore, $$R=X_{(n)}-X_{(1)}=\sum_{i=1}^n Y_i-Y_1=\sum_{i=2}^n Y_i$$
Hence, $$\mathbb E\left[R\right]=\sum_{i=2}^n \frac1{n-i+1}$$
And since $X_{(n)}=\sum\limits_{i=1}^n Y_i$, we also have $$\mathbb E\left[X_{(n)}\right]=\sum_{i=1}^n \mathbb E\left[Y_i\right]=\sum_{i=1}^n \frac1{n−i+1}=\sum_{i=1}^n \frac1{i}$$
Related threads:
Alternatively, we can proceed to find the expectation of $X_{(1)}$ and $X_{(n)}$ separately as you did. Clearly $X_{(1)}$ is exponential with mean $1/n$. And the density of $X_{(n)}$ is
$$f_{X_{(n)}}(x)=ne^{-x}(1-e^{-x})^{n-1}\mathbf1_{x>0}$$
For a direct calculation of the mean of $X_{(n)}$, we have
\begin{align}
\mathbb E\left[X_{(n)}\right]&=\int x f_{X_{(n)}}(x)\,dx
\\&=n\int_0^\infty xe^{-x}(1-e^{-x})^{n-1}\,dx
\\&=n\int_0^1(-\ln u)(1-u)^{n-1}\,du \tag{1}
\\&=n\int_0^1 -\ln(1-t)t^{n-1}\,dt \tag{2}
\\&=n\int_0^1 \sum_{j=1}^\infty \frac{t^j}{j}\cdot t^{n-1}\,dt \tag{3}
\\&=n\sum_{j=1}^\infty \frac1j \int_0^1 t^{n+j-1}\,dt \tag{4}
\\&=n\sum_{j=1}^\infty \frac1{j(n+j)}
\\&=\sum_{j=1}^\infty \left(\frac1j-\frac1{n+j}\right)
\\&=\sum_{j=1}^n \frac1j
\end{align}
$(1)$: Substitute $e^{-r}=u$.
$(2)$: Substitute $t=1-u$.
$(3)$: Use Maclaurin series expansion of $\ln(1-t)$ which is valid since $t\in (0,1)$.
$(4)$: Interchange integral and sum using Fubini/Tonelli's theorem.
We can also find the density of $R$ through the change of variables $(X_{(1)},X_{(n)})\mapsto (R,X_{(1)})$ and find $\mathbb E\left[R\right]$ directly by basically the same calculation as above.