Generating function for cubes of Harmonic numbers I

Let us define a following generating function: $${\bf H}^{(1,1,1)}_n(x) := \sum\limits_{m=1}^\infty [H_m]^3 \cdot \frac{x^m}{m^n}$$ Now, by using results from Generating function for cubes of Harmonic numbers we have found the following identity: \begin{eqnarray} &&{\bf H}^{(1,1,1)}_n(x) =\\ &&3! S_{n,3}(x) + 3! S_{n-1,4}(x) - 3 S_{n+1,2}(x)+Li_{n+3}(x)-3 S_{1,2}(x) Li_n(x)+\\ &&\frac{1}{2} \int\limits_0^x \frac{[Li_3(\xi)]^2}{\xi}\cdot \frac{[\log(x/\xi)]^{n-4}}{(n-4)!} d\xi+\\ &&3\int\limits_0^x \frac{Li_0(\xi) Li_1(\xi) Li_2(\xi)}{\xi}\cdot \frac{[\log(x/\xi)]^{n-1}}{(n-1)!} d\xi+\\ &&\frac{3}{2} \sum\limits_{l=1}^n \int\limits_0^x \frac{[Li_1(\xi)]^2}{\xi} Li_l(\xi) \cdot \frac{[\log(x/\xi)]^{n-l}}{(n-l)!} d\xi \end{eqnarray} where $S_{n,p}(x)$ is the Nielsen generalized poly-logarithm. The result above holds for $n\ge 4$ and $x\in(-1,1)$.

Now my question is how do we use this result in order to find closed form expressions for ${\bf H}^{(1,1,1)}_n(\pm1)$?

Let us look at the case $x=1$ in here. Note that the Nielsen generalized poly-logarithm at plus unity always reduces to single zeta values only. We know that from Compute an integral containing a product of powers of logarithms. for example.The only thing we need to do is to evaluate the three integrals. We do this in the usual way,i.e. by expanding the integrand in a multivariate series , integrating term by term and then using partial fraction decomposition to evaluate the sums. Now we state the results. We start from the easiest case and end with the hardest one. We have: \begin{eqnarray} &&\int\limits_0^1 \frac{[Li_3(\xi)]^2}{\xi} \cdot \frac{[\log(1/\xi)]^{n-4}}{(n-4)!} d\xi=\sum\limits_{\lambda_1\ge 1}\sum\limits_{\lambda_2\ge 1} \frac{1}{\lambda_1^3} \frac{1}{\lambda_2^3} \frac{1}{(\lambda_1+\lambda_2)^{n-3}}=\\ &&\sum\limits_{l_1=2}^3 \binom{n-1-l_1}{3-l_1} (-1)^{3-l_1} \zeta(l_1) \zeta(n-l_1+3) + \sum\limits_{l_1=2}^{n-3} \binom{n-1-l_1}{2} (-1)^3 \zeta(l_1,n-l_1+3) + \binom{n-2}{2} \left[ \zeta(n+2,1) + \zeta(n+3)\right] \end{eqnarray} The remaining two integrals require already a triple summation. We have: \begin{eqnarray} &&\int\limits_0^1 \frac{Li_0(\xi) Li_1(\xi) Li_2(\xi)}{\xi} \cdot \frac{[\log(1/\xi)]^{n-1}}{(n-1)!} d\xi= \sum\limits_{\lambda_1\ge 1}\sum\limits_{\lambda_2\ge 1}\sum\limits_{\lambda_3\ge 1} \frac{1}{\lambda_1^0} \frac{1}{\lambda_2^1} \frac{1}{\lambda_3^2} \frac{1}{(\lambda_1+\lambda_2+\lambda_3)^n}=\\ &&-\sum\limits_{i=1}^{n-1} \zeta(n+1-i,i,2)+\\ &&\zeta(2) \cdot\left[\zeta(n,1)+\zeta(n+1)\right]-\left[\zeta(2,n,1)+\zeta(2,n+1)\right]-\left[\zeta(n+2,1)+\zeta(n+3)\right] \end{eqnarray} and \begin{eqnarray} &&\int\limits_0^1 \frac{[Li_1(\xi)]^2}{\xi} \cdot Li_l(\xi) \cdot \frac{[\log(1/\xi)]^{n-l}}{(n-l)!} d\xi=\sum\limits_{\lambda_1\ge 1}\sum\limits_{\lambda_2\ge 1}\sum\limits_{\lambda_3\ge 1} \frac{1}{\lambda_1^1} \frac{1}{\lambda_2^1} \frac{1}{\lambda_3^l} \frac{1}{(\lambda_1+\lambda_2+\lambda_3)^{n-l+1}}=\\ &&-\sum\limits_{i=1}^{n-l} \sum\limits_{j=1}^l \zeta(n-l+2-i,i+j,l+1-j) \cdot\left(1_{j<l} + 2 \delta_{j,l}\right)+\\ &&\left.\sum\limits_{j=1}^l \left(\zeta(i+j+1,l+1-j)+\zeta(i+j,l+2-j) +\zeta(i+j,l+1-j,1)+\zeta(i+j,1,l+1-j)\right)\cdot\left(1_{j<l} + 2 \delta_{j,l}\right)\right|_{i=n-l+1} \end{eqnarray}