How to integrate $\int\frac{1}{\sqrt{1+x^3}}\mathrm dx$?

The first thing to do is to note that

$$x^3+1=(x+1)(x^2-x+1)$$

(one real and two complex conjugate roots). Using Jacobian elliptic functions requires having a quartic within the square root (the alternative of using Weierstrass elliptic functions is fine with square roots of cubics, but I'll leave that approach to someone else); the good thing is that by choosing a proper Möbius transformation, one can turn a cubic into a quartic (the algebraic geometers here might want to say a bit more than I have).

For the integral in question, the Möbius substitution needed is $x=\frac{-1+\sqrt{(-1)^2-(-1)+1}+(-1-\sqrt{(-1)^2-(-1)+1})v}{1+v}=\frac{2\sqrt{3}}{1+v}-(1+\sqrt{3})$; we then have

$$\int\frac{\mathrm dx}{\sqrt{x^3+1}}=-2\int\frac{\mathrm dv}{\sqrt{(1-v^2)(2\sqrt{3}-3+(2\sqrt{3}+3)v^2)}}$$

At this point, making use of the Jacobian elliptic function identity $\mathrm{sn}^2(u|m)+\mathrm{cn}^2(u|m)=1$ (nothing more than the usual Pythagorean identity in elliptic function garb), we could make either of the substitutions $v=\mathrm{sn}(u|m)$ or $v=\mathrm{cn}(u|m)$. The latter is a bit more convenient, since $\mathrm dv=-\mathrm{sn}(u|m)\mathrm{dn}(u|m)\mathrm du$, which can conveniently get rid of the negative sign in the integral. Thus, the integral turns into

$$2\int\frac{\mathrm{sn}(u|m)\mathrm{dn}(u|m)\mathrm du}{\sqrt{(1-\mathrm{cn}^2(u|m))(2\sqrt{3}-3+(2\sqrt{3}+3)\mathrm{cn}^2(u|m))}}$$

or (by using the Pythagorean identity)

$$2\int\frac{\mathrm{dn}(u|m)\mathrm du}{\sqrt{2\sqrt{3}-3+(2\sqrt{3}+3)\mathrm{cn}^2(u|m)}}$$

Here, one now chooses a proper value of $m$ such that the integrand reduces to a constant. Skipping the details, we let $m=\frac{2+\sqrt{3}}{4}$ such that

$$2\int\frac{\mathrm{dn}(u|m)\mathrm du}{\sqrt{2\sqrt{3}-3+(2\sqrt{3}+3)\mathrm{cn}^2(u|m)}}=\int\frac{\mathrm du}{\sqrt[4]{3}}$$

To undo the substitutions, we note that $u=F(\arccos(v)|m)$ and $v=\frac{2\sqrt{3}}{1+\sqrt{3}+x}-1$, giving the final result

$$\int\frac{\mathrm dx}{\sqrt{x^3+1}}=\frac1{\sqrt[4]{3}}F\left(\arccos\left(\frac{2\sqrt{3}}{1+\sqrt{3}+x}-1\right)\mid\frac{2+\sqrt{3}}{4}\right)+C$$

This result can be verified by differentiating the right hand side (remember that $\frac{\mathrm d}{\mathrm d\phi}F(\phi|m)=\frac1{\sqrt{1-m\sin^2\phi}}$) and noting that it is the same as the integrand.

For positive values of $x$, the integral is solved by $$ \int\sqrt{\dfrac1{1+x^3}}\mathrm \;dx = x \;_2F_1\left(\frac{1}{3},\frac{1}{2};\frac{4}{3},-x^3\right) $$ which contains a hypergeometric function. In general for positive $a$ $$ \int\sqrt{\dfrac1{a+x^3}}\mathrm \;dx = \frac{x}{\sqrt{a}} \;_2F_1\left(\frac{1}{3},\frac{1}{2};\frac{4}{3},-\frac{x^3}{a}\right) $$

Maybe it's off topic but I'd like to give the answer of Mathematica to the question "how to find the antiderivative in closed form, and how to prove it"

The antiderivative calculated by Mathematica is given by

ad = Integrate[1/Sqrt[1 + x^3 ], x]

(*
Out[347]= (2 (-1)^(1/6) Sqrt[-(-1)^(1/6) ((-1)^(2/3) + x)] Sqrt[
 1 + (-1)^(1/3) x + (-1)^(2/3) x^2]
  EllipticF[ArcSin[Sqrt[-(-1)^(5/6) (1 + x)]/3^(1/4)], (-1)^(1/3)])/(3^(
 1/4) Sqrt[1 + x^3])
*)

The result consists of an elliptic integral of the first kind and some square root factors. This looks a bit messy at first sight but it turns out that the square root factor cancel out.

Mathematica's (Full)Simplify does not simplify the expression, but my trick is to simplify the square of the expression in question.

FullSimplify[ad^2]

(*
Out[381]= (4 (-1)^(1/6)
  EllipticF[ArcSin[Sqrt[-(-1)^(5/6) (1 + x)]/3^(1/4)], (-1)^(1/3)]^2)/Sqrt[3]
*)

and take the square root afterwards (by hand). Then the antiderivative is

ad1 = (2 (-1)^(1/12)/3^(1/4)) EllipticF[
   ArcSin[Sqrt[-(-1)^(5/6) (1 + x)]/3^(1/4)], (-1)^(1/3)]

Mathematica's help tells us "EllipticF is the inverse of JacobiAmplitude. If $\phi=\mathrm{am}(u|m)$ then $u=F(\phi|m)$". Hence ad1 can can be expressed differently (as J. M. has ingeniously elaborated here).

Ok, now the proof, that ad1 is in fact an antiderivative (up to an additive constant, of course), differentiating ad1 with respect to x gives

D[ad, x]; (* -> horrible expression *)

Simplifying this, using my "square sandwiching trick" again, gives

FullSimplify[dad^2]

(*
Out[392]= 1/(1 + x^3)
*)

which completes the proof.