closed form solution for $\int\frac{I_{0}\left( x\right) }{e^{x}-1}dx$?

Question

Any one knows if there is a way to obtain a closed form solution to this integral?

$$ \int\frac{I_{0}\left( x\right) }{e^{x}-1}dx $$

where $I_{0}$ is modified Bessel function of first kind, also called BesselI.

The problem is when the denominator has the an extra term, such as -1 in it then no solution is found:

int(BesselI(0,x)/exp(x),x);

int(BesselI(0,x),x);

But

int(BesselI(0,x)/(exp(x)-1),x);

remains unevaluated.

Answer 1

$\int\dfrac{I_0(x)}{e^x-1}dx$

$=\int\sum\limits_{n=0}^\infty\dfrac{x^{2n}}{4^n(n!)^2(e^x-1)}dx$

$=\int\sum\limits_{n=0}^\infty\sum\limits_{k=0}^\infty\dfrac{B_kx^{2n+k-1}}{4^n(n!)^2k!}dx$ (with the formula in http://en.wikipedia.org/wiki/Bernoulli_number#Generating_function)

$=\int\left(\sum\limits_{n=0}^\infty\dfrac{x^{2n-1}}{4^n(n!)^2}+\sum\limits_{n=0}^\infty\sum\limits_{k=1}^\infty\dfrac{B_kx^{2n+k-1}}{4^n(n!)^2k!}\right)dx$

$=\int\left(\dfrac{1}{x}+\sum\limits_{n=1}^\infty\dfrac{x^{2n-1}}{4^n(n!)^2}+\sum\limits_{n=0}^\infty\sum\limits_{k=1}^\infty\dfrac{B_kx^{2n+k-1}}{4^n(n!)^2k!}\right)dx$

$=\ln x+\sum\limits_{n=1}^\infty\dfrac{x^{2n}}{2^{2n+1}(n!)^2n}+\sum\limits_{n=0}^\infty\sum\limits_{k=1}^\infty\dfrac{B_kx^{2n+k}}{4^n(n!)^2k!(2n+k)}+C$