Covariance of order statistics (uniform case)
I don't know if there is a smart way to solve the problem, but the standard way is the following:
For $1\le j< k\le n$, the joint density function should be $$ f_{X_{j},X_{_k}}=\frac{n!}{(j-1)!\,(k-j-1)!\,(n-k)!}\,x^{j-1}\,(y-x)^{k-1-j}\,(1-y)^{n-k} $$ where $0<x<y<1.$
Now to compute $E[X_{j}X_{_k}]=\int_0^1\mathbb{d}y\int_0^y xy\cdot f(x,y)\,\mathbb{d}x$, you need to be familiar with Beta function. Let $$C=\frac{n!}{(j-1)!\,(k-j-1)!\,(n-k)!},$$ so that $$\begin{eqnarray} E[X_{j}X_{_k}] &=&C\int_0^1\mathbb{d}y\int_0^y xy\cdot x^{j-1}(y-x)^{k-1-j}(1-y)^{n-k}\,\mathbb{d}x \\ &=&C\int_0^1\mathbb{d}y\int_0^y \left[\frac{x}{y}\right]^{j}\left[1-\frac{x}{y}\right]^{k-1-j}y^k(1-y^{n-k})\,\mathbb{d}x\\ &=&C\cdot B(j+1,k-j)\int_0^1y^{1+k}(1-y)^{n-k}\,\mathbb{d}y \\ &=&C \cdot B(j+1,k-j)\cdot B(k+2,n-k+1)\\ &=& \frac{j\cdot (k+1)}{(n+1)(n+2)}. \end{eqnarray}$$
Finally, $$ \begin{eqnarray}\operatorname{Cov}[X_{j}X_{_k}]&=&E[X_{j}X_{_k}]-E[X_{j}]\cdot E[X_{_k}]\\ &=& \frac{j\cdot (k+1)}{(n+1)(n+2)}-\frac{j\cdot k}{(n+1)^2} \\ &=& \frac{j\cdot (n+1-k)}{(n+1)^2(n+2)}. \end{eqnarray}$$
Thank you so much for posting this -- I too looked at this covariance integral (presented in David and Nagaraja's 2003 3rd edition Order Statistics text) and thought that it looked ugly. However, I fear that there may be a few small mistakes in your math on E(X_jX_k), assuming that I'm following you right. The joint density should have (j-1)! in the denominator instead of j! at the outset -- otherwise j! would entirely cancel out the j! in the numerator of Beta(j+1,k-j) instead of ending up with j in the numerator of the solution... Right?