Integrals of integer powers of dilogarithm function
There is a reasonably easy way to do these integrals. Above all, it is systematic. Each dilogarithm can be written as an iterated integral $$ L(x) = \int_0^x \frac{dt}{t}\int_0^t\frac{du}{1-u}. $$ Every power of a dilogarithm can then be expanded into a long product of integrals, giving, as an example: $$ \int_0^1 L(x)^2\,dx = \int_0^1dx \int_0^x \frac{dt_1}{t_1} \int_0^{t_1}\frac{du_1}{1-u_1} \int_0^x \frac{dt_2}{t_2} \int_0^{t_2}\frac{du_2}{1-u_2}. $$ The reason this is easy is that one integral representation of multiple zeta values has the form $$ \int_0^1 \frac{dw_1}{w_1-b_1}\int_0^{w_1}\frac{dw_2}{w_2-b_2}\cdots\int_0^{w_{n-1}}\frac{dw_n}{w_n-b_n}, $$ where $b_j\in\{0,1\}$.
In the case of $\int L(x)^n$, the domain of the integral is $$ 0< t_1,\ldots,t_n < x < 1, \qquad 0 < u_j < t_j, $$ which can be rearranged into a union of simplices of the form $$ 0<w_{2n}<\cdots<w_{1}<x<1, $$ (where $w_j$'s are all of $t_j$'s and $u_j$'s in some order), and the integral over each simplex has the form $$ \int_0^1dx\,\Omega, $$ where $\Omega$ is an iterated integral expression, starting with either $\int_0^x \frac{dw_1}{w_1}$ or $\int_0^x \frac{dw_1}{1-w_1}$.
It is easy to interchange the order of integration of $x$ and $w_1$, integrate over $x$, and be left with a very similar integral over another simplex with only the variables $w_2,\ldots,w_{2n}$ present.
Implementing this recursive procedure gives the following expressions for $I$: $$\begin{eqnarray} I_1 &=& \zeta(2) - 1 \\I_2 &=& -2\zeta(2) - 4\zeta(2,1) + 2\zeta(2,2) + 4\zeta(3,1) + 6 \\I_3 &=& 18\zeta(2) + 36\zeta(2,1) + 36\zeta(2,1,1) - 12\zeta(2,1,2) - 6\zeta(2,2) - 24\zeta(2,2,1) + 6\zeta(2,2,2) + 12\zeta(2,3,1) - 12\zeta(3,1) - 36\zeta(3,1,1) + 12\zeta(3,1,2) + 24\zeta(3,2,1) + 36\zeta(4,1,1) - 90 \\I_4 &=& -360\zeta(2) - 720\zeta(2,1) - 864\zeta(2,1,1) - 576\zeta(2,1,1,1) + 144\zeta(2,1,1,2) + 144\zeta(2,1,2) + 288\zeta(2,1,2,1) - 48\zeta(2,1,2,2) - 96\zeta(2,1,3,1) + 72\zeta(2,2) + 288\zeta(2,2,1) + 432\zeta(2,2,1,1) - 96\zeta(2,2,1,2) - 24\zeta(2,2,2) - 192\zeta(2,2,2,1) + 24\zeta(2,2,2,2) + 48\zeta(2,2,3,1) - 48\zeta(2,3,1) - 288\zeta(2,3,1,1) + 48\zeta(2,3,1,2) + 96\zeta(2,3,2,1) + 144\zeta(2,4,1,1) + 144\zeta(3,1) + 432\zeta(3,1,1) + 576\zeta(3,1,1,1) - 144\zeta(3,1,1,2) - 48\zeta(3,1,2) - 288\zeta(3,1,2,1) + 48\zeta(3,1,2,2) + 96\zeta(3,1,3,1) - 96\zeta(3,2,1) - 432\zeta(3,2,1,1) + 96\zeta(3,2,1,2) + 192\zeta(3,2,2,1) + 288\zeta(3,3,1,1) - 144\zeta(4,1,1) - 576\zeta(4,1,1,1) + 144\zeta(4,1,1,2) + 288\zeta(4,1,2,1) + 432\zeta(4,2,1,1) + 576\zeta(5,1,1,1) + 2520 \end{eqnarray}$$
I computed these values up to $n=8$, but the expressions are long ($I_5$ has 428 terms), and the integers involved become prohibitively large, so I didn't include them here.
For small $n$, at least up to $n=4$, I found numerically that the above expressions can be written in terms of ordinary zeta values, as follows: $$\begin{eqnarray}\def\tfrac#1#2{{\textstyle\frac{#1}{#2}}} I_1 &=& \zeta(2) - 1 \\I_2 &=& -4 \zeta (3)-2 \zeta (2)+\tfrac{5}{2} \zeta (4)+6 \\I_3 &=& -12 \zeta (3) \zeta (2)+6 \zeta (5)+36 \zeta (3)+18 \zeta (2)+\tfrac{57}{2} \zeta (4)+\tfrac{35}{8} \zeta (6)-90 \\I_4 &=& 24 \zeta (5) \zeta (2)+144 \zeta (3) \zeta (2)+84 \zeta (3) \zeta (4)-195 \zeta (7)-648 \zeta (5)+144 \zeta (3)^2-720 \zeta (3)-360 \zeta (2)-774 \zeta (4)-\frac{11}{2} \zeta (6)+\frac{175}{24} \zeta (8)+2520 \end{eqnarray}$$ I checked $n\leq3$ using known identities satisfied by MZVs, such as those here. For $I_4$, while this is not a proof that the expression is correct, this entire procedure reduces the problem of evaluating the integrals $I_n$ to the problem of applying known identities to multiple zeta values, which is more straightforward.
A quick integer relation search for a closed form for $I_5$ gave me nothing.
There may or may not be a closed form for $I_5$ (I believe this is not a settled question), but there is no systematic way to reduce multiple zeta values to ordinary zeta values. For example, weight 10, which is what $I_5$ requires, is the first weight at which known relations between depth-two zeta values are insufficient to reduce them to ordinary zeta values ("Euler Sums and Contour Integral Representations" by Flajolet and Salvy). Because here the depth can be up to five, this problem is harder.
Update (by editor after 6 years): Settled now. By further computation one have
$$\small\int_0^1 \text{Li}_2(x){}^5 \, dx=1440 \zeta(6,2)+240 \zeta (3)^3+120 \pi ^2 \zeta (3)^2-4320 \zeta (3)^2-\frac{23 \pi ^6 \zeta (3)}{189}-110 \pi ^4 \zeta (3)-600 \pi ^2 \zeta (3)-9360 \zeta (5) \zeta (3)+25200 \zeta (3)-\frac{127 \pi ^4 \zeta (5)}{3}-540 \pi ^2 \zeta (5)+30600 \zeta (5)-\frac{325 \pi ^2 \zeta (7)}{2}+14445 \zeta (7)+\frac{16495 \zeta (9)}{3}+\frac{\pi ^{10}}{7776}+\frac{51137 \pi ^8}{45360}+\frac{1907 \pi ^6}{126}+335 \pi ^4+2100 \pi ^2-113400$$
Here $ \zeta(6,2)$ is irreducible. In fact all $I_n$ can be reduced to MZVs. The generalized algorithm I used also gives the proof to closed-form of $I_4$ above.
Allow me to present another approach. Though lengthier and more tedious, this method is independent of results derived from previous answers. The main idea is to expand ${\rm Li}_2(x)$ as a series and integrate term by term.
For $n=1$, \begin{align} \mathcal{I}_1 &=\int^1_0{\rm Li}_2(x) \ {\rm d}x\\ &=\int^1_0\sum^\infty_{j=1}\frac{x^j}{j^2}{\rm d}x\\ &=\sum^\infty_{j=1}\frac{1}{j^2(j+1)}\\ &=\sum^\infty_{j=1}\frac{1}{j^2}-\frac{1}{j}+\frac{1}{j+1}\\ \end{align} We see that \begin{align} \sum^\infty_{k=1}\frac{1}{k^2} &=\frac{\pi^2}{6}\\ \sum^\infty_{k=1}\left(\frac{1}{k+1}-\frac{1}{k}\right) &=\sum^\infty_{k=2}\frac{1}{k}-\sum^\infty_{k=1}\frac{1}{k}\\ &=-1 \end{align} Therefore $$\mathcal{I}_1=\frac{\pi^2}{6}-1$$
For $n=2$, \begin{align} \mathcal{I}_2 &=\int^1_0{\rm Li}^2_2(x) \ {\rm d}x\\ &=\int^1_0\sum^\infty_{j=1}\frac{x^j}{j^2}\sum^\infty_{k=1}\frac{x^k}{k^2}{\rm d}x\\ &=\sum^\infty_{j=1}\frac{1}{j^2}\sum^\infty_{k=1}\frac{1}{k^2(k+j+1)}\\ &=\sum^\infty_{j=1}\frac{1}{j^2(j+1)}\sum^\infty_{k=1}\frac{1}{k^2}-\frac{1}{k(j+1)}+\frac{1}{(k+j+1)(j+1)}\\ &=\frac{\pi^2}{6}\left(\frac{\pi^2}{6}-1\right)+\sum^\infty_{j=1}\frac{1}{j^2(j+1)^2}\sum^\infty_{k=1}\frac{1}{k+j+1}-\frac{1}{k}\\ &=\frac{\pi^4}{36}-\frac{\pi^2}{6}+\color{red}{\sum^\infty_{j=1}\frac{2H_{j+1}}{j}}\color{green}{-\sum^\infty_{j=1}\frac{2H_{j+1}}{j+1}}\color{blue}{-\sum^\infty_{j=1}\frac{H_{j+1}}{j^2}}\color{purple}{-\sum^\infty_{j=1}\frac{H_{j+1}}{(j+1)^2}}\\ \end{align} It is trivial to show that \begin{align} \sum^\infty_{j=1}\frac{H_j}{j^2} &=2\zeta(3) \end{align} We see that \begin{align} \color{red}{\sum^\infty_{j=1}\frac{2H_{j+1}}{j}} &=\sum^\infty_{j=1}\frac{2H_j}{j}+\frac{2}{j(j+1)}\\ &=\color{red}{\sum^\infty_{j=1}\frac{2H_j}{j}+2}\\ \color{green}{-\sum^\infty_{j=1}\frac{2H_{j+1}}{j+1}} &=\color{green}{-\sum^\infty_{j=1}\frac{2H_j}{j}+2}\\ \color{blue}{-\sum^\infty_{j=1}\frac{H_{j+1}}{j^2}} &=-\sum^\infty_{j=1}\frac{H_{j}}{j^2}-\sum^\infty_{j=1}\frac{1}{j^2(j+1)}\\ &=\color{blue}{-2\zeta(3)+1-\frac{\pi^2}{6}}\\ -\color{purple}{\sum^\infty_{j=1}\frac{H_{j+1}}{(j+1)^2}} &=-\sum^\infty_{j=1}\frac{H_{j}}{j^2}+1\\ &=\color{\purple}{-2\zeta(3)+1} \end{align} Therefore $$\mathcal{I}_2=\frac{\pi^4}{36}-4\zeta(3)-\frac{\pi^2}{3}+6$$
For $n=3$, \begin{align} \mathcal{I}_3 &=\int^1_0{\rm Li}^3_2(x) \ {\rm d}x\\ &=\int^1_0\sum^\infty_{j=1}\frac{x^j}{j^2}\sum^\infty_{k=1}\frac{x^k}{k^2}\sum^\infty_{m=1}\frac{x^m}{m^2}{\rm d}x\\ &=\sum^\infty_{j=1}\frac{1}{j^2}\sum^\infty_{k=1}\frac{1}{k^2}\sum^\infty_{m=1}\frac{1}{m^2(m+j+k+1)}\\ &=\sum^\infty_{j=1}\frac{1}{j^2}\sum^\infty_{k=1}\frac{1}{k^2(j+k+1)}\sum^\infty_{m=1}\frac{1}{m^2}-\frac{1}{m(j+k+1)}+\frac{1}{(m+j+k+1)(j+k+1)}\\ &=\frac{\pi^2}{6}\mathcal{I}_2-\sum^\infty_{j=1}\frac{1}{j^2}\sum^\infty_{k=1}\frac{H_{k+j+1}}{k^2(k+j+1)^2}\\ &=\frac{\pi^2}{6}\mathcal{I}_2-\sum^\infty_{j=1}\frac{1}{j^2(j+1)^2}\sum^\infty_{k=1}\left(\frac{H_{k+j+1}}{k^2}+\frac{H_{k+j+1}}{(k+j+1)^2}\right)+2\sum^\infty_{j=1}\frac{1}{j^2(j+1)^3}\sum^\infty_{k=1}\left(\frac{H_{k+j+1}}{k}-\frac{H_{k+j+1}}{k+j+1}\right) \end{align} The first inner sum is \begin{align} \sum^\infty_{k=1}\frac{H_{k+j+1}}{k^2} &=2\zeta(3)+\sum^\infty_{k=1}\sum^{j+1}_{m=1}\frac{1}{k^2(k+m)}\\ &=2\zeta(3)+\sum^\infty_{k=1}\sum^{j+1}_{m=1}\frac{1}{mk^2}-\frac{1}{m^2k}+\frac{1}{m^2(k+m)}\\ &=2\zeta(3)+\frac{\pi^2}{6}H_{j+1}-\sum^{j+1}_{m=1}\frac{1}{m^2}\sum^\infty_{k=1}\left(\frac{1}{k}-\frac{1}{k+m}\right)\\ &=2\zeta(3)+\frac{\pi^2}{6}H_{j+1}-\sum^{j+1}_{k=1}\frac{H_k}{k^2}\\ \end{align} The second sum is \begin{align} \sum^\infty_{k=1}\frac{H_{k+j+1}}{(k+j+1)^2} &=2\zeta(3)-\sum^{j+1}_{k=1}\frac{H_k}{k^2} \end{align} Adding them up gives $$4\zeta(3)+\frac{\pi^2}{6}H_{j+1}-2\sum^{j+1}_{k=1}\frac{H_k}{k^2}$$ The third sum is \begin{align} \sum^\infty_{k=1}\frac{H_{k+j+1}}{k} &=\sum^\infty_{k=1}\frac{H_k}{k}+\sum^\infty_{k=1}\sum^{j+1}_{m=1}\frac{1}{k(k+m)}\\ &=\sum^\infty_{k=1}\frac{H_{k}}{k}+\sum^\infty_{k=1}\sum^{j+1}_{m=1}\frac{1}{mk}-\frac{1}{m(k+m)}\\ &=\sum^\infty_{k=1}\frac{H_{k}}{k}+\sum^{j+1}_{k=1}\frac{H_k}{k}\\ \end{align} The fourth sum is \begin{align} \sum^\infty_{k=1}\frac{H_{k+j+1}}{k+j+1} &=\sum^\infty_{k=1}\frac{H_{k}}{k}-\sum^{j+1}_{k=1}\frac{H_k}{k} \end{align} Subtracting them gives $$H_{j+1}^2+H_{j+1}^{(2)}$$ So $$\mathcal{I}_3=\frac{\pi^2}{6}\mathcal{I}_2-\sum^\infty_{j=1}\frac{1}{j^2(j+1)^2}\left(4\zeta(3)+\frac{\pi^2}{6}H_{j+1}-2\sum^{j+1}_{k=1}\frac{H_k}{k^2}\right)+2\sum^\infty_{k=1}\frac{H_{j+1}^2+H_{j+1}^{(2)}}{j^2(j+1)^3}$$ I will try to complete the evaluation of $\mathcal{I}_3$ over the next few days, and add in the computation of $\mathcal{I}_4$ using this method if I have the time. It is clear that $\mathcal{I}_3$ can be reduced to a bunch of Euler sums, which are all rather manageable to compute.
Integrate by parts twice. \begin{align} \int^1_0{\rm Li}_2^3(x){\rm d}x &=\left[x{\rm Li}^3_2(x)\right]^1_0+\int^1_03{\rm Li}^2_2(x)\ln(1-x){\rm d}x\\ &=\frac{\pi^6}{216}-\frac{\pi^4}{12}+6\color\red{\int^1_0\frac{\left[(x-1)\ln(1-x)-x\right]{\rm Li}_2(x)\ln(1-x)}{x}{\rm d}x}\\ &=\frac{\pi^6}{216}-\frac{\pi^4}{12}+6\zeta(5)+\frac{2\pi^4}{5}+36\zeta(3)+3\pi^2-2\pi^2\zeta(3)-90\\ &=\frac{\pi^6}{216}+\frac{19\pi^4}{60}+6\zeta(5)+36\zeta(3)+3\pi^2-2\pi^2\zeta(3)-90\\ \end{align} The red integral was evaluated in your previous answer.