Consider the hypergeometric equation with parameters $(a,b,c)=\left(\frac16,\frac12,\frac13\right)$, and build from its two canonical solutions near $z=0$ the vector
$$\vec{y}(z)=\left(\begin{array}{c}
y_1 \\ y_2
\end{array}\right)=\left(\begin{array}{c}
_2F_1(a,b;c;z) \\ z^{1-c}{}_2F_1(a-c+1,b-c+1;2-c;z)
\end{array}\right).\tag{1}$$
This is a single-valued vector function on $\mathbb{C}\backslash\{(-\infty,0]\cup[1,\infty)\}$. Its analytic continuation along a closed loop $\gamma$ gives rise to monodromy representation of $\pi_1(\mathbb{C}\backslash\{0,1\})$:
$$ \gamma\mapsto M_{[\gamma]},\qquad y(\gamma z)=M_{[\gamma]}y(z).$$
The monodromy group $G\subset GL(2,\mathbb{C})$ of the hypergeometric equation is generated by two matrices corresponding to simple loops around $0$ and $1$. In the case we are interested in these matrices are explicitly given by
$$M_0=\left(\begin{array}{cc} 1 & 0 \\ 0 & e^{-2\pi i /3}\end{array}\right),\qquad
M_1=C\left(\begin{array}{cc} 1 & 0 \\ 0 & e^{2\pi i /3}\end{array}\right)C^{-1},\tag{2}$$
where the connection matrix $C=\left(\begin{array}{cc} 1 & 2^{\frac43} \\ -2^{\frac83} & 8\end{array}\right)$. If $G$ is finite, then $\vec{y}(z)$ has a finite number of branches, and moreover (Schwarz 1872), is algebraic.
It is not difficult to check that the monodromy group $G$ generated by $M_0$, $M_1$ from (2) is indeed finite. In particular, note that
$$M_0^3=M_1^3=I,\qquad M_1^2=-M_0M_1M_0, $$ $$M_1M_0M_1=-M_0^2,\qquad M_1M_0^2M_1=M_0^2M_1M_0^2.$$
It turns out that $G$ has order $24$ and is isomorphic to the binary tetrahedral group:
$$G\cong 2T=\langle s,t\,|\,(st)^2=s^3=t^3\rangle, $$
where the generators can be identified as $s=M_0M_1M_0M_1M_0$, $t=M_1M_0M_1M_0M_1$.
Corollary: The hypergeometric functions in (1) are algebraic.
Algebraic solutions of the hypergeometric equations are classified by the so-called Schwarz table, and have been studied by many mathematicians, see e.g. the bibliography in this paper by R.Vidunas. Their explicit construction is somewhat involved but relatively straightforward - at least when the corresponding algebraic curve has genus $0$ (the genus can be determined independently from the Riemann-Hurwitz formula).
In our case the task simplifies even more as our parameter values can be obtained from the genus $0$ tetrahedral formula (2.4) of the above mentioned paper by a combination of a linear trasformation (sending $\frac56$ to $\frac43-\frac56=\frac12$) and differentiation (transforming $\frac43$ into $\frac13$). The result is
$$_2F_1\left(\frac16,\frac12;\frac13;-\frac{r(r+2)^3}{(r+1)(1-r)^3}\right)=\frac{\sqrt{1-r^2}}{2r+1}.$$
Corollary: The antiderivative $\displaystyle\int \mathcal{R}\left(x,y(x)\right)dx$, where $y(x)={}_2F_1\left(\frac16,\frac12;\frac13;-x\right)$ and $\mathcal{R}(x,y)$ is rational in both arguments, can be expressed in terms of elementary functions.
Example:
The transformation $r\mapsto x(r)=\frac{r(r+2)^3}{(r+1)(1-r)^3}$ bijectively maps $(0,1)$ to $(0,\infty)$, and therefore the initial integral becomes
\begin{align}\mathcal{I}&=\int_0^1 \left(\frac{\sqrt{1-r^2}}{2r+1}\right)^{12}\left(\frac{r(r+2)^3}{(r+1)(1-r)^3}\right)'dr=\\&=
2\int_0^1 \frac{(1+r)^4(1-r)^2(r+2)^2}{(2r+1)^{10}} dr=\\&=\frac{80\,663}{153\,090}.
\end{align}
Edited for a more concise derivation:
Define $\mathcal{I}$ to be the value of the definite integral,
$$\mathcal{I}:=\int_{0}^{1}\frac{\ln^{2}{\left(x\right)}}{\sqrt{x^{2}-x+1}}\,\mathrm{d}x.\tag{1}$$
The definite integral $\mathcal{I}$ is found to have as an approximate numerical value
$$\mathcal{I}\approx2.10029.\tag{2}$$
We begin by transforming the integral via an Euler substitution of the first kind:
$$\sqrt{x^{2}-x+1}=x+t\implies x=\frac{1-t^{2}}{1+2t},\tag{3}$$
$$\implies\mathrm{d}x=\frac{\left(-2\right)\left(1+t+t^{2}\right)}{\left(1+2t\right)^{2}}\,\mathrm{d}t,$$
$$\implies\sqrt{x^{2}-x+1}=\frac{1-t^{2}}{1+2t}+t=\frac{1+t+t^{2}}{1+2t}.$$
Under the transformation $(3)$, the integral $\mathcal{I}$ becomes
$$\begin{align}
\mathcal{I}
&=\int_{0}^{1}\frac{\ln^{2}{\left(x\right)}}{\sqrt{x^{2}-x+1}}\,\mathrm{d}x\\
&=2\int_{0}^{1}\frac{\ln^{2}{\left(\frac{1-t^{2}}{1+2t}\right)}}{1+2t}\,\mathrm{d}t;~~~\small{\left[\sqrt{x^{2}-x+1}=x+t\right]}\\
&=2\int_{0}^{1}\frac{\left[\ln{\left(1-t\right)}+\ln{\left(1+t\right)}-\ln{\left(1+2t\right)}\right]^{2}}{1+2t}\,\mathrm{d}t.\tag{4}\\
\end{align}$$
Next, using the algebraic identity
$$\left(a+b-c\right)^{2}=a^{2}+2b^{2}-\left(a-b\right)^{2}+\left(a-c\right)^{2}-2bc,$$
with $a=\ln{\left(1-t\right)}\land b=\ln{\left(1+t\right)}\land c=\ln{\left(1+2t\right)}$, we may expand the integral $\mathcal{I}$ as a sum of five simpler logarithmic integrals:
$$\begin{align}
\mathcal{I}
&=2\int_{0}^{1}\frac{\left[\ln{\left(1-t\right)}+\ln{\left(1+t\right)}-\ln{\left(1+2t\right)}\right]^{2}}{1+2t}\,\mathrm{d}t\\
&=2\int_{0}^{1}\frac{\ln^{2}{\left(1-t\right)}}{1+2t}\,\mathrm{d}t+4\int_{0}^{1}\frac{\ln^{2}{\left(1+t\right)}}{1+2t}\,\mathrm{d}t\\
&~~~~~-2\int_{0}^{1}\frac{\left[\ln{\left(1-t\right)}-\ln{\left(1+t\right)}\right]^{2}}{1+2t}\,\mathrm{d}t\\
&~~~~~+2\int_{0}^{1}\frac{\left[\ln{\left(1-t\right)}-\ln{\left(1+2t\right)}\right]^{2}}{1+2t}\,\mathrm{d}t\\
&~~~~~-4\int_{0}^{1}\frac{\ln{\left(1+t\right)}\ln{\left(1+2t\right)}}{1+2t}\,\mathrm{d}t\\
&=2\int_{0}^{1}\frac{\ln^{2}{\left(1-t\right)}}{1+2t}\,\mathrm{d}t+4\int_{0}^{1}\frac{\ln^{2}{\left(1+t\right)}}{1+2t}\,\mathrm{d}t\\
&~~~~~-2\int_{0}^{1}\frac{\ln^{2}{\left(\frac{1-t}{1+t}\right)}}{1+2t}\,\mathrm{d}t+2\int_{0}^{1}\frac{\ln^{2}{\left(\frac{1-t}{1+2t}\right)}}{1+2t}\,\mathrm{d}t\\
&~~~~~-4\int_{0}^{1}\frac{\ln{\left(1+t\right)}\ln{\left(1+2t\right)}}{1+2t}\,\mathrm{d}t.\tag{5}\\
\end{align}$$
Now, we'll find it convenient to introduce the following two auxiliary functions:
$$J{\left(z\right)}:=\int_{0}^{1}\frac{\ln^{2}{\left(y\right)}}{1+zy}\,\mathrm{d}y;~~~\small{z\ge-1},\tag{6a}$$
and
$$H{\left(a,c\right)}:=\int_{0}^{1}\frac{\ln^{2}{\left(1+ay\right)}}{1+cy}\,\mathrm{d}y;~~~\small{a\ge-1\land c>-1}.\tag{6b}$$
As we shall see, the integral $\mathcal{I}$ can be expressed entirely in terms of the auxiliary functions $J$ and $H$, and hence, entirely in terms of elementary functions and the standard polylogarithms:
$$\begin{align}
\mathcal{I}
&=2\int_{0}^{1}\frac{\ln^{2}{\left(1-t\right)}}{1+2t}\,\mathrm{d}t+4\int_{0}^{1}\frac{\ln^{2}{\left(1+t\right)}}{1+2t}\,\mathrm{d}t\\
&~~~~~-2\int_{0}^{1}\frac{\ln^{2}{\left(\frac{1-t}{1+t}\right)}}{1+2t}\,\mathrm{d}t+2\int_{0}^{1}\frac{\ln^{2}{\left(\frac{1-t}{1+2t}\right)}}{1+2t}\,\mathrm{d}t\\
&~~~~~-\int_{0}^{1}\frac{4\ln{\left(1+t\right)}\ln{\left(1+2t\right)}}{1+2t}\,\mathrm{d}t\\
&=2\int_{0}^{1}\frac{\ln^{2}{\left(u\right)}}{3-2u}\,\mathrm{d}u;~~~\small{\left[1-t=u\right]}\\
&~~~~~+4\,H{\left(1,2\right)}\\
&~~~~~-\int_{0}^{1}\frac{4\ln^{2}{\left(v\right)}}{\left(3-v\right)\left(1+v\right)}\,\mathrm{d}v;~~~\small{\left[\frac{1-t}{1+t}=v\right]}\\
&~~~~~+2\int_{0}^{1}\frac{\ln^{2}{\left(w\right)}}{1+2w}\,\mathrm{d}w;~~~\small{\left[\frac{1-t}{1+2t}=w\right]}\\
&~~~~~\small{-\left(\left[\ln{\left(1+t\right)}\ln^{2}{\left(1+2t\right)}\right]_{0}^{1}-\int_{0}^{1}\frac{\ln^{2}{\left(1+2t\right)}}{1+t}\,\mathrm{d}t\right)};~~~\small{I.B.P.s}\\
&=\frac23\int_{0}^{1}\frac{\ln^{2}{\left(u\right)}}{1-\frac23u}\,\mathrm{d}u+4\,H{\left(1,2\right)}\\
&~~~~~-\frac13\int_{0}^{1}\frac{\ln^{2}{\left(v\right)}}{1-\frac13v}\,\mathrm{d}v-\int_{0}^{1}\frac{\ln^{2}{\left(v\right)}}{1+v}\,\mathrm{d}v+2\,J{\left(2\right)}\\
&~~~~~-\ln{(2)}\ln^{2}{(3)}+H{\left(2,1\right)}\\
&=\frac23\,J{\left(-\frac23\right)}+4\,H{\left(1,2\right)}-\frac13\,J{\left(-\frac13\right)}\\
&~~~~~-J{\left(1\right)}+2\,J{\left(2\right)}-\ln{(2)}\ln^{2}{(3)}+H{\left(2,1\right)}.\tag{7}\\
\end{align}$$
The function $H$ can be expressed in terms of $J$, dilogarithms, and elementary functions. Define $\gamma:=\frac{a-c}{c}$. For $0<a\land0<c$, we have $-1<\gamma$ and
$$\begin{align}
H{\left(a,c\right)}
&=\int_{0}^{1}\frac{\ln^{2}{\left(1+ay\right)}}{1+cy}\,\mathrm{d}y\\
&=\int_{\frac{1}{1+a}}^{1}\frac{\ln^{2}{\left(x\right)}}{x\left[c+\left(a-c\right)x\right]}\,\mathrm{d}x;~~~\small{\left[\frac{1}{1+ay}=x\right]}\\
&=\frac{1}{c}\int_{\frac{1}{1+a}}^{1}\frac{\ln^{2}{\left(x\right)}}{x}\,\mathrm{d}x-\frac{a-c}{c}\int_{\frac{1}{1+a}}^{1}\frac{\ln^{2}{\left(x\right)}}{c+\left(a-c\right)x}\,\mathrm{d}x\\
&=\frac{1}{3c}\ln^{3}{\left(1+a\right)}-\frac{a-c}{c^{2}}\int_{\frac{1}{1+a}}^{1}\frac{\ln^{2}{\left(x\right)}}{1+\left(\frac{a-c}{c}\right)x}\,\mathrm{d}x\\
&=\frac{1}{3c}\ln^{3}{\left(1+a\right)}-\frac{\gamma}{c}\int_{\frac{1}{1+a}}^{1}\frac{\ln^{2}{\left(x\right)}}{1+\gamma x}\,\mathrm{d}x\\
&=\frac{1}{3c}\ln^{3}{\left(1+a\right)}-\frac{\gamma}{c}\int_{0}^{1}\frac{\ln^{2}{\left(x\right)}}{1+\gamma x}\,\mathrm{d}x\\
&~~~~~+\frac{\gamma}{c}\int_{0}^{\frac{1}{1+a}}\frac{\ln^{2}{\left(x\right)}}{1+\gamma x}\,\mathrm{d}x\\
&=\frac{1}{3c}\ln^{3}{\left(1+a\right)}-\frac{\gamma}{c}J{\left(\gamma\right)}\\
&~~~~~+\frac{\gamma}{c\left(1+a\right)}\int_{0}^{1}\frac{\ln^{2}{\left(\frac{w}{1+a}\right)}}{1+\left(\frac{\gamma}{1+a}\right)w}\,\mathrm{d}w;~~~\small{\left[\left(1+a\right)x=w\right]}\\
&=\frac{1}{3c}\ln^{3}{\left(1+a\right)}-\left(\frac{a-c}{c^{2}}\right)J{\left(\frac{a-c}{c}\right)}\\
&~~~~~\small{+\frac{\gamma}{c\left(1+a\right)}\int_{0}^{1}\frac{\ln^{2}{\left(w\right)}-2\ln{\left(w\right)}\ln{\left(1+a\right)}+\ln^{2}{\left(1+a\right)}}{1+\left(\frac{\gamma}{1+a}\right)w}\,\mathrm{d}w}\\
&=\frac{1}{3c}\ln^{3}{\left(1+a\right)}-\left(\frac{a-c}{c^{2}}\right)J{\left(\frac{a-c}{c}\right)}\\
&~~~~~+\frac{\gamma}{c\left(1+a\right)}\int_{0}^{1}\frac{\ln^{2}{\left(w\right)}}{1+\left(\frac{\gamma}{1+a}\right)w}\,\mathrm{d}w\\
&~~~~~-\frac{2}{c}\ln{\left(1+a\right)}\int_{0}^{1}\frac{\left(\frac{\gamma}{1+a}\right)\ln{\left(w\right)}}{1+\left(\frac{\gamma}{1+a}\right)w}\,\mathrm{d}w\\
&~~~~~+\frac{1}{c}\ln^{2}{\left(1+a\right)}\int_{0}^{1}\frac{\left(\frac{\gamma}{1+a}\right)}{1+\left(\frac{\gamma}{1+a}\right)w}\,\mathrm{d}w\\
&=\frac{1}{3c}\ln^{3}{\left(1+a\right)}-\left(\frac{a-c}{c^{2}}\right)J{\left(\frac{a-c}{c}\right)}\\
&~~~~~+\frac{\gamma}{c\left(1+a\right)}\,J{\left(\frac{\gamma}{1+a}\right)}\\
&~~~~~+\frac{2}{c}\ln{\left(1+a\right)}\int_{0}^{1}\frac{\ln{\left(1+\left(\frac{\gamma}{1+a}\right)w\right)}}{w}\,\mathrm{d}w\\
&~~~~~+\frac{1}{c}\ln^{2}{\left(1+a\right)}\ln{\left(1+\left(\frac{\gamma}{1+a}\right)\right)}\\
&=\frac{1}{3c}\ln^{3}{\left(1+a\right)}-\left(\frac{a-c}{c^{2}}\right)J{\left(\frac{a-c}{c}\right)}\\
&~~~~~+\frac{a-c}{c^{2}\left(1+a\right)}\,J{\left(\frac{a-c}{c\left(1+a\right)}\right)}\\
&~~~~~-\frac{2}{c}\ln{\left(1+a\right)}\operatorname{Li}_{2}{\left(-\left(\frac{\gamma}{1+a}\right)\right)}\\
&~~~~~+\frac{1}{c}\ln^{2}{\left(1+a\right)}\ln{\left(\frac{a\left(1+c\right)}{\left(1+a\right)c}\right)}\\
&=\frac{a-c}{c^{2}\left(1+a\right)}\,J{\left(\frac{a-c}{c\left(1+a\right)}\right)}-\frac{a-c}{c^{2}}\,J{\left(\frac{a-c}{c}\right)}\\
&~~~~~-\frac{2}{c}\ln{\left(1+a\right)}\operatorname{Li}_{2}{\left(\frac{c-a}{c\left(1+a\right)}\right)}\\
&~~~~~+\frac{1}{3c}\ln^{3}{\left(1+a\right)}+\frac{1}{c}\ln^{2}{\left(1+a\right)}\ln{\left(\frac{a\left(1+c\right)}{\left(1+a\right)c}\right)}.\tag{8}\\
\end{align}$$
The function $J$ reduces to the trilogarithm. For $-1\le z\land z\neq0$,
$$\begin{align}
J{\left(z\right)}
&=\int_{0}^{1}\frac{\ln^{2}{\left(y\right)}}{1+zy}\,\mathrm{d}y\\
&=-\frac{2}{z}\operatorname{Li}_{3}{\left(-z\right)}.\tag{9}\\
\end{align}$$
Thus, continuing from where we left off at the last line of $(7)$,
$$\begin{align}
\mathcal{I}
&=\frac23\,J{\left(-\frac23\right)}-\frac13\,J{\left(-\frac13\right)}-J{\left(1\right)}+2\,J{\left(2\right)}\\
&~~~~~+4\,H{\left(1,2\right)}+H{\left(2,1\right)}-\ln{(2)}\ln^{2}{(3)}\\
&=2\operatorname{Li}_{3}{\left(\frac23\right)}-2\operatorname{Li}_{3}{\left(\frac13\right)}+2\operatorname{Li}_{3}{\left(-1\right)}-2\operatorname{Li}_{3}{\left(-2\right)}\\
&~~~~~\small{+4\left[-\operatorname{Li}_{3}{\left(\frac14\right)}+\operatorname{Li}_{3}{\left(\frac12\right)}-\ln{(2)}\operatorname{Li}_{2}{\left(\frac14\right)}-\frac56\ln^{3}{(2)}+\frac12\ln^{2}{(2)}\ln{(3)}\right]}\\
&~~~~~\small{+\left[2\operatorname{Li}_{3}{\left(-1\right)}-2\operatorname{Li}_{3}{\left(-\frac13\right)}-\ln{(9)}\operatorname{Li}_{2}{\left(-\frac13\right)}-\frac23\ln^{3}{(3)}+\ln{(4)}\ln^{2}{(3)}\right]}\\
&~~~~~-\ln{(2)}\ln^{2}{(3)}\\
&=2\operatorname{Li}_{3}{\left(\frac23\right)}-2\operatorname{Li}_{3}{\left(\frac13\right)}-2\operatorname{Li}_{3}{\left(-\frac13\right)}-2\operatorname{Li}_{3}{\left(-2\right)}\\
&~~~~~-4\operatorname{Li}_{3}{\left(\frac14\right)}+4\operatorname{Li}_{3}{\left(\frac12\right)}+4\operatorname{Li}_{3}{\left(-1\right)}\\
&~~~~~-4\ln{(2)}\operatorname{Li}_{2}{\left(\frac14\right)}-2\ln{(3)}\operatorname{Li}_{2}{\left(-\frac13\right)}\\
&~~~~~-\frac{10}{3}\ln^{3}{(2)}+2\ln^{2}{(2)}\ln{(3)}+\ln{(2)}\ln^{2}{(3)}-\frac23\ln^{3}{(3)}.\blacksquare\\
\end{align}$$
For $z\notin[1,\infty)$,
$$\small{\operatorname{Li}_{3}{\left(z\right)}=-\operatorname{Li}_{3}{\left(\frac{z}{z-1}\right)}-\operatorname{Li}_{3}{\left(1-z\right)}+\frac{\ln^3{\left(1-z\right)}}{6}-\frac{\ln{\left(z\right)}\ln^2{\left(1-z\right)}}{2}+\frac{\pi^2}{6}\ln{\left(1-z\right)}+\zeta{(3)}}.$$
For $0<z<1$,
$$\small{\operatorname{Li}_{3}{\left(z\right)}+\operatorname{Li}_{3}{\left(1-z\right)}+\operatorname{Li}_{3}{\left(\frac{z}{z-1}\right)}=\frac{\ln^3{\left(1-z\right)}}{6}-\frac{\ln{\left(z\right)}\ln^2{\left(1-z\right)}}{2}+\frac{\pi^2}{6}\ln{\left(1-z\right)}+\zeta{(3)}}.$$
Setting $z=\frac13$,
$$\small{\operatorname{Li}_{3}{\left(\frac13\right)}+\operatorname{Li}_{3}{\left(\frac23\right)}+\operatorname{Li}_{3}{\left(-\frac12\right)}=\frac{\ln^3{\left(\frac23\right)}}{6}-\frac{\ln{\left(\frac13\right)}\ln^2{\left(\frac23\right)}}{2}+\zeta{(2)}\ln{\left(\frac23\right)}+\zeta{(3)}}.$$
$$\small{\operatorname{Li}_{3}{\left(\frac13\right)}+\operatorname{Li}_{3}{\left(\frac23\right)}+\operatorname{Li}_{3}{\left(-\frac12\right)}=\frac{\ln^{3}{(2)}-3\ln{(2)}\ln^{2}{(3)}+2\ln^{3}{(3)}}{6}+\zeta{(2)}\ln{\left(\frac23\right)}+\zeta{(3)}}$$
$$\small{\operatorname{Li}_{3}{\left(\frac13\right)}+\operatorname{Li}_{3}{\left(\frac23\right)}+\operatorname{Li}_{3}{\left(-\frac12\right)}=\frac{\ln^{2}{\left(\frac32\right)}\ln{\left(18\right)}}{6}-\zeta{(2)}\ln{\left(\frac32\right)}+\zeta{(3)}}$$
Best Answer
What about the Laplace transform? By using it, we have that our integral equals:
$$ I= \frac{1}{36}\int_{0}^{+\infty}\left(1-3 s+6 s^2-6 s^3 \psi'(1+s)\right)\,ds$$ and in this form Mathematica is perfectly able to state that $I=\frac{\zeta(3)}{\color{red}{8}\pi^2}$.
I just used: $$\mathcal{L}^{-1}\left(\frac{1}{x^4}\right)=\frac{s^3}{6},\qquad \mathcal{L}\left(1-\frac{x}{2}+\frac{x^2}{12}-\frac{x}{e^x-1}\right)=\frac{1-3 s+6 s^2}{6 s^3}-\psi'(1+s) $$ together with: $$ \int_{0}^{+\infty}f(x)g(x)\,dx = \int_{0}^{+\infty}(\mathcal{L} f)(s)(\mathcal{L}^{-1}g)(s)\,ds.$$