Ordered statistics: find $E(Y(i) Y(j))$ Let $Y1 ,…., Yn$ be independent, uniformly distributed random variables on the interval [0,1]. 
Find: $Var[Y_{(j)}- Y_{(i)}]$, the variance of the difference between two ordered statistics 
Where i and j are integers with 1≤i < j ≤n.
My Steps are:
1) find the joint density function. My answer:
$$\frac{n!}{(i-1)!(j-i-1)!(n-j)!}(y_{(i)})^{i-1}(y_{(j)}-y_{(i)})^{j-1-i}(1-y_{(j)})^{n-j}$$ 
where $0< y_{(i)}< y_{(j)}<1$
2) Find $E[Y_{(i)}Y_{(j)}]$ from first principle
$$a\int_{0}^{1}\int_{y(_{j})}^{1}y_{(j)}(y_{(j)}-y_{(i)})^{(j-i-1)}(1-y_{(j)})^{n-j}y_{i}^{i} dy_{(i)}dy_{(j)}$$
where $$a = \frac{n!}{(i-1)!(j-i-1)!(n-j)!}$$
This is where I am stuck. I am not sure how to integrate this integrand. Any tips/hints would be great. My aim is to convert this to a Beta form.
3)  Find $\operatorname{Cov}[Y_{(i)}, Y_{(j)}]$ using $\operatorname{Cov}[Y_{(i)},Y_{(j)}] = E[Y_{(i)}Y_{(j)}] - E[Y_{(i)}]E[Y_{(j)}]$. But this depends on step 2.
4)  Find the variance using the generic $\operatorname{Var}[A-B]=\operatorname{Var}[A]+\operatorname{Var}[B]-2\operatorname{Cov}[A,B]$ . However, this all depends on part 3 and the steps before it.
Please provide some guidance as to whether my overall approach is correct and what can I do to overcome step 2. 
Thanks. 
Edit: sorry for step 2 -> it should be $E[Y_{(i)}Y_{(j)}]$ not $E[Y_{(i)}-Y_{(j)}]$. Typing error.
 A: Why don't you try and find $E[Y_{(j)}-Y_{(i)}]$ by using the linearity of expectation.  You don't need the double integral.
$$E[Y_{(j)}-Y_{(i)}]=E[Y_{(j)}]-E[Y_{(i)}]$$
Thus, the pdf for the $k$th order statistic is
$$f_k(y)=\frac{n!}{(k-1)!(n-k)!}y^{k-1}(1-y)^{n-k}dy$$
$$E(Y_k)=\int_0^1{y\frac{n!}{(k-1)!(n-k)!}y^{k-1}(1-y)^{n-k}}dy$$
$$=\int_0^1{\frac{n!}{(k-1)!(n-k)!}y^{k}(1-y)^{n-k}}dy$$
From here it is a bunch of integration by parts, or you rewrite the integral in one step:
$$=\frac{k}{n+1}\int_0^1{\frac{(n+1)!}{(k)!(n-k)!}y^{k}(1-y)^{n-k}}dy$$
Now the integral expression is just the beta distribution with parameters $\alpha=k+1$and $\beta=n-r+1$.  Thus the entire integral is just 1!  Hence
$$E[Y_{(j)}-Y_{(i)}]=E[Y_{(j)}]-E[Y_{(i)}]=\frac{j}{n+1}-\frac{i}{n+1}=\frac{j-i}{n+1}$$
Try seeing how far you can get now that you have this part....
A: Here's the solution generated by mathStatica/Mathematica …
Parent pdf:  let $Y$ ~ Uniform(0,1) with pdf $f(y)$, which we enter as:
f = 1;        domain[f] = {y, 0, 1};

The joint pdf of the (i, j)th order statistics, say pdf $g(y_{i}, y_{j})$, is given by:
g = OrderStat[{i, j}, f]


with domain statement:
  domain[g] = OrderStatDomain[{i, j}, f];

Part 3: $\operatorname{Cov}[Y_{(i)}, Y_{(j)}]$:

Part 4: The variance you seek is:

All done.
