What is wrong with this approach to find Expected value of distance between $X, Y \in$ Uniform $(0, 1)$? Let $X$ and $Y$ be i.i.d. random variables with Uniform $(0, 1)$ continuous distribution.
The problem is to find the expected value of the distance between X and Y.
My reasoning was, for all $(x, y) \in X \times Y$, any distance must be $0 \le |x - y| \le 1$, so we can consider two cases, $0 \le y - x \le 1$ and $0 \le x - y \le 1$. 
So to find the expectated value of the distance, I would take the sum of the two integrals,
$$ \int_0^1 \int_{y-1}^1 (y - x) \; dx \; dy $$
$$ \int_0^1 \int_{x-1}^1 (x- y) \; dy \; dx $$
The first and second integrals each give $1/3$, so the answer I reached is $2/3$, whereas the correct answer for expected distance is $1/3$ and uses integrals over different limits of integration. I know what the correct answer says, but where am I going wrong here?
 A: The setting up of the integrals was not correct. If $X$ and $Y$ are our random variables, we want to find 
$$\iint_{S} |x-y|\,dxdy,$$
where $S$ is the $1\times 1$ square. 
Draw that square. Now it is useful to break up this integral into two parts, by drawing the line $y=x$, which slices our square into two parts: (i) where $y\ge x$ and (ii) where $y\le x$.
To do (i), express as an iterated integral. If you wish to integrate first with respect to $x$, note that $x$ travels from $0$ to $y$. The travels of $x$ can be seen clearly in the picture. So for (i) we need
$$\int_{y=0}^1\left(\int_{x=0}^y(y-x)\,dx    \right)\,dy.$$
Because of symmetry, (ii) will give the same result. 
There are a couple of problems with your integrals. First note that there has not even been an integration with respect to $y$. Secondly, the setting up of the limits has no clear connection to the geometry. 
A: Edit: my previous answer was actually not correct, so I removed it completely.
I had said that you need to divide your integral results by 2, since they each represented half of the domain space. This is not correct; if the integrals' bounds are set up correctly and we are actually integrating over the correct areas, the fact that they individually account for only half of the domain will naturally be accounted-for.
The only error in your original approach is the setup of the integrals. So, draw the image of $[0,1]\times[0,1]$ as a visual aid, and work from there. The integrals corresponding to your two cases should be $$\int_{0}^{1}\int_{0}^{y}(y-x)dxdy$$ for the upper triangle, where $y>x$, and $$\int_{0}^{1}\int_{0}^{x}(x-y)dydx$$, where $x>y$. By symmetry these integrals will have the same value, so double the calculation from each. (each has value 1/6).
Sorry for the confusion.
A: This is not really an answer to the question "Where am I going wrong here?" 
but rather an alternative way of approaching the problem to be solved that 
gives the answer quite easily.
$Z = |X-Y|$ is a non-negative random variable taking on values in $[0,1]$, and hence
$$E[Z] = \int_0^{\infty} P\{Z > z\}\,\mathrm dz = \int_0^1 P\{Z > z\}\,\mathrm dz.$$
But, for $0 \leq z \leq 1$, $P\{Z > z\}$ is the probability that the random point
$(X,Y)$ lies in one of two right-triangular regions of area $\frac{1}{2}(1-z)^2$ each.
This is easiest to figure out if you follow André Nicolas's suggestion to draw
a sketch: the triangles have vertices $(z,0), (1,0), (1,1-z)$ and $(0,z), (0,1),(1-z,1)$
respectively. Hence,
$$E[Z] = \int_0^1 P\{Z > z\}\,\mathrm dz = \int_0^1(1-z)^2\,\mathrm dz
= \left. \frac{-1}{3}(1-z)^3\right|_0^1 = \frac{1}{3}.$$
