Using Harmonic polylogarithms to solve a partial differential equation

Question

I am trying to solve the following partial differential equation:

Eq1:=$f''[x] - ((2*(2*x^2 - 1))/(x*(x^2 - 1)))*f'[x] + 4/(x^2 - 1) f[x] = HPL[{4}, x]/(x*(1 - x)*(1 + x))$

It is told to us that mathematica cannot solve this equation with DSolve[]. We need to use the Harmonic polylogarithms (HPL). We needed to istall this package: https://www.physik.uzh.ch/data/HPL/, https://arxiv.org/abs/hep-ph/0507152.

Now, whee it comes to solving the PDE, I use the following steps:

1. I find the solution to the hom. PDE.
2. Once I have all the solution:$f_{h,i}(x),i=1,...m$ for the particular solution I can say:

$$f_p(x)=\sum_{i=1}^n f_{h,i}(x)\int \frac{g(y)W_i(y)}{W(y)}$$, where

$W(x)=det(\frac{\partial^k}{\partial x^k}f_{h,m}(x))$, k=0,...n-1 (for a n-th order PDE) and m=1,....n

and

$W_i(x)=(-1)^{i+n}det(\frac{\partial^k}{\partial x^k}f_{h,m}(x))$.

Initially I solve the hom. PDE by using DSolve[Eq1].

Normally, Mathematica can integrate such functions only in their standard form, i.e:

Integrate[ HPL[{0,1,-1},x]/(1-x) , {x,0,xx}]

For more complex ones, it is needed to simplify the expression, or rather, the division:

HPL[n_,xx_]/(1+xx_) -> Integrate[ HPL[n,xx]/(1+xx) , {xx,0,x}]

But, I cannot understand how, can I ultimately have in my integral the W(..) terms, as these resemble the result of a matrix.

How can I find the particular solution in this case?