How to evaluate $\int_0^1\frac{\log^2(1+x)}x\mathrm dx$?
The similar integral
$$ \int_0^1\frac{\log^2(1-x)}x\mathrm dx=2\zeta(3) $$
is evaluated in this blog post using the substitution $u=-\log(1-x)$:
$$ \begin{align} \int_0^1\frac{\log^2(1-x)}x\mathrm dx &= \int_0^\infty\frac{u^2}{1-\mathrm e^{-u}}\mathrm e^{-u}\,\mathrm du \\ &= \int_0^\infty u^2\sum_{n=1}^\infty\mathrm e^{-nu}\mathrm du \\ &= \sum_{n=1}^\infty\int_0^\infty u^2\mathrm e^{-nu}\mathrm du \\ &= \sum_{n=1}^\infty\frac2{n^3} \\ &= 2\zeta(3)\;. \end{align} $$
Analogously substituting $u=\log(1+x)$ in the present integral leads to an integral up to $\log2$ that can be evaluated in terms of polylogarithms evaluated at $\frac12$:
$$ \begin{align} &\int_0^{\log2}\frac{\log^2(1+x)}x\mathrm dx \\ =& \int_0^{\log2}\frac{u^2}{\mathrm e^u-1}\mathrm e^u\,\mathrm du \\ =& \int_0^{\log2}\frac{u^2}{1-\mathrm e^{-u}}\mathrm du \\ =& \int_0^{\log2} u^2\sum_{n=0}^\infty\mathrm e^{-nu}\mathrm du \\ =& \sum_{n=0}^\infty\int_0^{\log2} u^2\mathrm e^{-nu}\mathrm du \\ =& \sum_{n=0}^\infty\int_0^{\log2} u^2\mathrm e^{-nu}\mathrm du \\ =& \frac13\log^32+\sum_{n=1}^\infty\frac1n\left(-2^{-n}\log^22+2\int_0^{\log2} u\mathrm e^{-nu}\mathrm du\right) \\ =& \frac13\log^32+\sum_{n=1}^\infty\left(-\frac1n2^{-n}\log^22+\frac2{n^2}\left(-2^{-n}\log2+\int_0^{\log2}\mathrm e^{-nu}\mathrm du\right)\right) \\ =& \frac13\log^32+\sum_{n=1}^\infty\left(-\frac1n2^{-n}\log^22-\frac2{n^2}2^{-n}\log2-\frac2{n^3}\left(2^{-n}-1\right)\right) \\ =& \def\Li{\operatorname{Li}} \frac13\log^32-\Li_1\left(\frac12\right)\log^22-2\Li_2\left(\frac12\right)\log2-2\Li_3\left(\frac12\right)+2\zeta(3) \\ =& \frac13\log^32-\log2\log^22-2\left(\frac{\pi^2}{12}-\frac12\log^22\right)\log2-2\left(\frac16\log^32-\frac{\pi^2}{12}\log2+\frac78\zeta(3)\right)+2\zeta(3) \\ =& \frac{\zeta(3)}4\;. \end{align} $$
Not only is this a rather complicated derivation of a much simpler result; it also looks as if the polylogarithm values may have been obtained using the present integral in the first place.
We can combine the present integral with the similar integral
$$ \int_0^1\frac{\log^2(1-x)}x\mathrm dx=2\zeta(3) $$
(see my other answer) into
$$ \begin{align} \int_0^1\frac{\log^2(1+x)}x\mathrm dx-\int_0^1\frac{\log^2(1-x)}x\mathrm dx &= \int_0^1\frac{\log^2(1+x)}x\mathrm dx+\int_{-1}^0\frac{\log^2(1+x)}x\mathrm dx \\ &= \int_{-1}^1\frac{\log^2(1+x)}x\mathrm dx\;. \end{align} $$
Then we can complete the contour of integration by a semicircle in the upper half complex plane:
$$ \begin{align} \int_{-1}^1\frac{\log^2(1+x)}x\mathrm dx &= \oint\frac{\log^2(1+z)}z\mathrm dz-\int\frac{\log^2\left(1+\def\e{\mathrm e^{\mathrm i\phi}}\e\right)}{\e}\mathrm d\e \\ &= -\mathrm i\int_0^\pi\log^2\left(1+\e\right)\mathrm d\phi\;, \end{align} $$
where the integral over the closed contour vanishes since there are no poles inside the contour.
We know that the imaginary part of this expression vanishes, since it sums to zero with a real integral, so we only have to evaluate the real part:
$$ \begin{align} -\mathrm i\int_0^\pi\log^2\left(1+\e\right)\mathrm d\phi &= \Re\left(-\mathrm i\int_0^\pi\log^2\left(1+\e\right)\mathrm d\phi\right) \\ &= \Re\left(-\mathrm i\int_0^\pi\left(\log\left|1+\e\right|+\mathrm i\arg\left(1+\e\right)\right)^2\mathrm d\phi\right) \\ &= 2\int_0^\pi\log\left|1+\e\right|\arg\left(1+\e\right)\mathrm d\phi \\ &= 2\int_0^\pi\frac12\log\left(2+2\cos\phi\right)\frac\phi2\mathrm d\phi \\ &= \frac12\int_0^\pi\log\left(\left(\mathrm e^{\mathrm i\phi/2}+\mathrm e^{-\mathrm i\phi/2}\right)^2\right)\phi\mathrm d\phi \\ &= \int_0^\pi\log\left(\mathrm e^{\mathrm i\phi/2}+\mathrm e^{-\mathrm i\phi/2}\right)\phi\,\mathrm d\phi \\ &= \Re\int_0^\pi\left(\frac{\mathrm i\phi}2+\log\left(1+\mathrm e^{-\mathrm i\phi}\right)\right)\phi\,\mathrm d\phi \\ &= \Re\int_0^\pi\sum_{n=1}^\infty\frac{(-1)^{n+1}\mathrm e^{-\mathrm in\phi}}n\phi\,\mathrm d\phi \\ &= \sum_{n=1}^\infty\frac{-1+(-1)^n}{n^3} \\ &= -\zeta(3)-\eta(3) \\ &= -\zeta(3)-\frac34\zeta(3) \\ &= -\frac74\zeta(3)\;. \end{align} $$
The desired integral is the sum of the two results:
$$ \begin{align} \int_0^1\frac{\log^2(1+x)}x\mathrm dx &= \int_0^1\frac{\log^2(1-x)}x\mathrm dx-\mathrm i\int_0^\pi\log^2\left(1+\e\right)\mathrm d\phi \\ &= 2\zeta(3)-\frac74\zeta(3) \\ &= \frac{\zeta(3)}4 \;. \end{align} $$
This raises the question whether there's a deeper reason for both of these seemingly quite different integrals to evaluate to a multiple of $\zeta(3)$.
Let,
$\displaystyle A=\int_0^1 \dfrac{\ln\left(1+x\right)^2}{x}dx$
$\displaystyle B=\int_0^1 \dfrac{\ln\left(1-x\right)^2}{x}dx$
Perform the change of variable $y=1-x$,
$\displaystyle B=\int_0^1 \dfrac{\ln\left(x\right)^2}{1-x}dx$
Using Taylor expansion,
$\displaystyle B=2\zeta(3)$
$\displaystyle C=\int_0^1 \dfrac{\ln\left(\dfrac{1-x}{1+x}\right)^2}{x}dx$
Perform the change of variable $y=\dfrac{1-x}{1+x}$,
$\displaystyle C=2\int_0^1 \dfrac{\ln\left(x\right)^2}{1-x^2}dx$
Using Taylor expansion,
$C=\dfrac{7}{2}\zeta(3)$
$\displaystyle D=\int_0^1 \dfrac{\ln\left(1-x^2\right)^2}{x}dx$
Perform the change of variable $y=x^2$,
$\begin{align}\displaystyle D&=\dfrac{1}{2}\int_0^1 \dfrac{\ln\left(1-x\right)^2}{x}dx\\ &=\dfrac{1}{2}B\\ &=\zeta(3) \end{align}$
Since $(a+b)^2+(a-b)^2=2a^2+2b^2$ therefore,
$C+D=2A+2B$
Therefore,
$\begin{align}A&=\dfrac{1}{2}\left(C+D-2B\right)\\ &=\dfrac{1}{2}\left(\dfrac{7}{2}\zeta(3)+\zeta(3)-4\zeta(3)\right)\\ &=\boxed{\dfrac{1}{4}\zeta(3)} \end{align}$
$\newcommand{\+}{^{\dagger}} \newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle} \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack} \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,} \newcommand{\dd}{{\rm d}} \newcommand{\down}{\downarrow} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,{\rm e}^{#1}\,} \newcommand{\fermi}{\,{\rm f}} \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{{\rm i}} \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow} \newcommand{\isdiv}{\,\left.\right\vert\,} \newcommand{\ket}[1]{\left\vert #1\right\rangle} \newcommand{\ol}[1]{\overline{#1}} \newcommand{\pars}[1]{\left(\, #1 \,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}} \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,} \newcommand{\sech}{\,{\rm sech}} \newcommand{\sgn}{\,{\rm sgn}} \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert} \newcommand{\wt}[1]{\widetilde{#1}}$ \begin{align}&\color{#c00000}{\int_{0}^{1}{\ln^{2}\pars{1 + x} \over x}\,\dd x} =\int_{1}^{2}{\ln^{2}\pars{x} \over x - 1}\,\dd x =\int_{1}^{1/2}{\ln^{2}\pars{1/x} \over 1/x - 1}\,\pars{-\,{\dd x \over x^{2}}} \\[3mm]&=\int_{1/2}^{1}{\ln^{2}\pars{x} \over x\pars{1 - x}}\,\dd x =\int_{1/2}^{1}{\ln^{2}\pars{x} \over x}\,\dd x + \int_{1/2}^{1}{\ln^{2}\pars{x} \over 1 - x}\,\dd x \\[3mm]&={1 \over 3}\,\ln^{3}\pars{2} +\color{#66f}{\sum_{n = 0}^{\infty}\int_{1/2}^{1}\ln^{2}\pars{x}x^{n}\,\dd x} \end{align}
\begin{align}&\color{#66f}{\sum_{n = 0}^{\infty}\int_{1/2}^{1}\ln^{2}\pars{x}x^{n} \,\dd x} =\left.\partiald[2]{}{\mu}\sum_{n = 1}^{\infty}\int_{1/2}^{1}x^{\mu - 1} \,\dd x\,\right\vert_{\,\mu\ =\ n} =\left.\partiald[2]{}{\mu}\sum_{n = 1}^{\infty} {1 - 2^{-\mu} \over \mu}\,\right\vert_{\,\mu\ =\ n} \\[3mm]&=2\sum_{n = 1}^{\infty}{1 \over n^{3}} -2\sum_{n = 1}^{\infty}{\pars{1/2}^{n} \over n^{3}} -2\ln\pars{2}\sum_{n = 1}^{\infty}{\pars{1/2}^{n} \over n^{2}} -\ln^{2}\pars{2}\sum_{n = 1}^{\infty}{\pars{1/2}^{n} \over n} \\[3mm]&=2\zeta\pars{3} - 2{\rm Li}_{3}\pars{\half} -2\ln\pars{2}{\rm Li}_{2}\pars{\half} -\ln^{2}\pars{2}{\rm Li}_{1}\pars{\half} \end{align}
From this link \begin{align} {\rm Li}_{1}\pars{\half} &= \ln\pars{2} \\[1mm] {\rm Li}_{2}\pars{\half} &= {\pi^{2} \over 12} - \half\,\ln^{2}\pars{2} \\[1mm] {\rm Li}_{3}\pars{\half} &= {1 \over 6}\,\ln^{3}\pars{2} -{\pi^{2} \over 12}\,\ln\pars{2} + {7 \over 8}\,\zeta\pars{3} \end{align}