This sum is a function of the form $f(x_1,\ldots,x_n)=x^TMx+\lambda^Tx$ with $x=[x_1,\ldots,x_n]^T$ and vector $\lambda\in\mathbb{R}^n$. If you carry out the calculations it takes the form
$$f(x)=\big(\sum_{i=1}^n{x_i}\big)^2-n\sum_{i=1}^n{x_i^2}-\sum_{i=1}^n{(n-2i+1)x_i}.\qquad \qquad (1)$$
The key result is that we can write (1) equivalently as
$$f(x)=-n\sum_{i=1}^n{\Big(x_i-\frac{1}{n}\sum_{j=1}^{n}{x_j}+\frac{n-2i+1}{2n}\Big)^2}+\frac{1}{4n}{\sum_{i=1}^n(n-2i+1)^2} \qquad\qquad (2)$$
which achieves maximum whenever
$$x_i=\frac{1}{n}\sum_{j=1}^{n}{x_j}-\frac{n-2i+1}{2n}\qquad \qquad\qquad\qquad\qquad\qquad\qquad(3)$$
This maximum value is attainable if $0\leq x_i\leq 1$ or equivalently if
$\frac{n-2i+1}{2n}\leq \frac{1}{n}\sum_{j=1}^{n}{x_j}\leq\frac{n-2i+1}{2n}+ 1$ for all $i=1,\ldots,n$
i.e. if
$$\frac{n-1}{2n}\leq \frac{1}{n}\sum_{j=1}^{n}{x_j}\leq\frac{n+1}{2n}$$
Selecting $\frac{1}{n}\sum_{j=1}^{n}{x_j}=\frac{n+1}{2n}$ we obtain from (3) $x_i=i/n$. Now if we replace (3) in (2) the maximum value of
$$\frac{1}{4n}{\sum_{i=1}^n(n-2i+1)^2}=\frac{n^2-1}{12}$$
is derived.
An interesting observation on the problem is that there is a infinite number of points that achieve this maximum value even in $[0,1]^n$. This is a consequence of (3) where one can see that for each value of the sum $\sum_{i=1}^n{x_i}$ in the interval $[(n-1)/2,(n+1)/2]$ different values of the optimal $x_i$ are obtained. They only have to differ by $1/n$.