The integral I need to evaluate is:
$$ \int_x^{\infty} \frac{t^n}{e^{t} -1} dt $$
After some research I found a paper saying,
The numerical values of the two integrals [...] are easily calculated either directly or through series expansion.
It's easy to evaluate, right? So I try plugging this into my 32bit 7.7.0.471 (R2008b) MATLAB program and it cannot integrate it. Fortunately they actually cited this and I was lead to a mathematics table which stated, $$\int_x^{\infty} \frac{t^n dt}{e^t -1} = \sum_{k=1}^{\infty} \exp^{-kx} ( \frac{x^n}{k} + \frac{n x^{n-1}}{k^2} + \frac{ (n)(n-1) x^{n-2}}{k^3} + \cdots +\frac{n!}{k^{n+1}} )$$
I'm using $n=2$ so my actual expression is $$\int_x^{\infty} \frac{t^2 dt}{e^t -1} = \sum_{k=1}^{\infty} \exp(-kx) ( \frac{x^2}{k} + \frac{2 x^{1}}{k^2} + \frac{ (2)(1) x^{0}}{k^3})$$
The problem is, however, that my $x$ is very large for the environment, for example $10^{19}$, and the precision for $\exp$ causes it to be evaluated as 0. If I try to scale this on a $\log$ basis then I can't seem to get it converge even after several thousand summations. My software is also too old to be compatible with newer contributions which would solve my problem by allowing me to arbitrarily increase the precision.
This paper is from the mid 90's, why did they say it is so easy? Is there a better way, mathematically, to attack this that a novice like myself can understand and implement? How does one even find this series? It seems crazy.
I did find a paper that seems to be doing what I want, although I am not sure how to implement it.