You can discretize the equation by approximating the derivatives by finite differences and use numerical integration. Once you do that, you'll have a nonlinear system of equations that may be numerically solved by some fixed point algorithm (for instance Newton's method).
Setting up a grid of equally spaced points in $[0,1]$, say $u_i = i h$, with $h = \frac 1n$, the equation can be approximated by (I'm switching $\xi$ to $y$):
$$
-6 \dfrac{y_{i+1}-2y_1+y_{i-1}}{h^2} \sum_{i=0}^n \frac{w_i}{y_i^2}+4m\left( -2 \dfrac{y_{i+1}-2y_1+y_{i-1}}{h^2} + 2a^2c^2\sum_{i=0}^n \frac{w_i}{y_i^3}\right) = \frac{6}{y_i^3}\left(\frac{w_0(y_1-y_0)^2}{h^2} + \sum_{i=1}^{n-1} w_i \left(\frac{y_{i+1}-y_{i-1}}{2h}\right)^2 + \frac{w_n(y_n-y_{n-1})^2}{h^2}\right)
$$
The coefficients $w_i$ depend on the quadrature rule that you use to approximate the integrals. In the case of the trapezoidal rule, they are given by $\frac h2, h, h, \cdots, h, \frac h2$. These equations are for the interior nodes. For the boundary nodes you must add two additional equations: One of them must be $y_0 - y_n = 0$, enforcing the periodic boundary conditions; The other can be a modification os the equation for interior nodes, but using forward (if using $u_0$) or backward (if using $u_n$) finite differences, in stead of the centred ones used before.
Finally, you just need to use your favorite solver to solve the system of nonlinear equations for the unknowns $y_0, \cdots, y_n$. Since you mentioned Mathematica, I suggest $\texttt{FindRoot[]}$.
Below you can find the numerical solution using 4 grid nodes and 40 grid nodes. All the constants were set to 1 and, just out of laziness, I've set the last equation to $y_n=1$, instead of changing to backward differences.

