I did it with Lagrange multipliers. The problem is equivalent to maximizing $f(x_1,...,x_n) = \displaystyle\sum_{i<j} x_i x_j(x_i^2 + x_j^2) = \sum_i (1-x_i)x_i^3$, with the constraint that $\sum x_i = 1$.
Using the Lagrangian, at a maximum that isn't on the boundary $ 4x_i^3 - 3x_i^2 = 4x_j^3 - 3x_j^2\hspace {10pt} \forall i,j$.
$f(\dfrac 1 n, ..., \dfrac 1 n) = \dfrac {n-1} {n^3} < f(\dfrac 1 2,\dfrac 1 2, 0, ..., 0)$ for $n \geq 3$.
The function $g(x) = 4x^3 - 3x^2$ has a local maximum at 0 and a local minimum at $x=\dfrac 1 2$.
So the only other possibility for a maximum is a point of the type $(w,z,...,z)$ where $w > \dfrac 1 2, 1-w = (n-1)z$.
Let $n_1 = n-1$.
In that case, $f(w,z,...,z) = n_1(1-z)z^3 + n_1z(1-n_1z)^3$.
But $n_1(1-z)z^3 < (1-n_1z)(n_1z)^3$ for $n_1 \geq 2$, since $z < \dfrac 1 {2n_1}$.
So $f(w,z,...,z) < f(w,1-w)$.
But for $n=2$, at the maximum $4((1-w)^2+(1-w)w + w^2) = 3((1-w)+w)= 3$, which is only true if $w=\dfrac 1 2$.
So the max possible value of $f$ for any $n$ is $\dfrac 1 8$, at a point of the type $(\dfrac 1 2, \dfrac 1 2, 0, ..., 0)$.