Continued fraction expansion related to exponential generating function A recent SciComp.SE Question motivates us to ask for a nice continued fraction expansion of the following Maclaurin series:
$$ f(x) = \sum_{n=0}^\infty \frac{B_n\; x^{n+3}}{n! (n+3)} = \int_0^x \frac{t^3}{e^t - 1} dt $$
where $B_n$ are the (first) Bernoulli numbers.
Although the positive power of $t$ in the numerator of the integrand removes the singularity at the origin, the exponential function has complex period $2\pi i$.  Therefore the radius of convergence for the above series is $2\pi$, the distance from the origin to the nearest singularities.
Where power series convergence is always limited by the proximity of singularities, continued fraction expansions typically have larger domains of convergence, perhaps including the whole real line.  Such a circumstance weighs in favor of rational over polynomial approximations over a widened real interval, and nice continued fractions often give apt starting values for a minimax rational approximation.
 A: Background:
It is evident from the integral form that the function $f(x)$ will have a triple root at $x=0$, that it is monotone increasing, and (due to the integrand's denominator growing exponentially) finitely bounded above (by exactly $\pi^4/15$ as it turns out).
Here are a few fairly precise function values that illustrate this trend:
$$\begin{array}{|r|c|} 
\hline \text{x} & \text{f(x) (shown rounded)} \\ \hline
1.0 & 0.22480518802594 \\ \hline
2.0 & 1.17634259660700 \\ \hline
3.0 & 2.55221845329080 \\ \hline
4.0 & 3.87705416153119 \\ \hline
5.0 & 4.89989215833058 \\ \hline
6.0 & 5.58585538083094 \\ \hline
7.0 & 6.00316896121307 \\ \hline
8.0 & 6.23962379489192 \\ \hline
9.0 & 6.36657389887547 \\ \hline
10.0 & 6.43192189678183 \\ \hline
\end{array}$$
It was also evident in from the Question itself that the function has a power series expansion about the origin, but that the radius of convergence is only $2\pi$ owing to the nearest singularities (poles).  Since these singularities "accumulate" at infinity in the complex plane, it cannot be expected that a convergent power series representation can be centered "at infinity", i.e. a Laurent series.
Nevertheless it may prove desirable to investigate that in more detail, after first exploring the possibility of converting the Maclaurin series into a continued fraction.  Below the third-order Debye function being considered is graphed alongside its truncated Maclaurin series polynomials of order 7, 15, and 23 for visual consistency:

Fig. 1 Debye function of 3rd order and its Maclaurin polynomials
Notice the diminishing returns near the threshold radius of convergence, that the improvement in approximation going from order 15 to order 23 is slight in comparison with the earlier improvement from order 7 to order 15.  As nonconstant polynomials inevitably veer off to plus or minus infinity, it is natural to hope for better convergence from rational expressions (which allow for the horizontal asymptote seen here).
Graphic produced using gnuplot interface from Maxima; I'll post the Maxima code as a GitHub gist once I've polished it...
A: Your integral term $$f(t)= t^3/(e^t-1)$$ can be expressed as the continued fraction $$f(t) = t^2/(1+t/(2+t/( -{3 \over 2}+t/(-8+t/({5 \over 6}+t/(18+t/(-{7 \over 12} + ... $$
$$f(t) ={ t^2 \over1+ K_{i=2}^\infty{ t \over c_i} }\ where \ c_i =  \begin{cases}i\ even \ c_i=(-1)^{(i+2) \over 2}2({i \over 2})^2 \\ i \ odd \ \ \,c_i = (-1)^{(i-1) \over 2}{i \over { {i-1} \over 2}{{i+1} \over 2} }
 \end{cases} \ $$
Wallis' Fundamental Recurrence Formulas provide an easy way to calculate your optimal rational polynomials.
$$ f(t) \approx {h_i \over k_i}$$
$$ h_i = 1 ,\ 0, \ t^2, \ \ \ \   ( 2t^2 ), \ \  (t^3-3t^2), \ \  (-6t^3 +24t^2), 
 (t^4 -8t^2+20t^2),\ ...$$
$$ k_i = 0, \ 1, \ 1, \  \ (t+2), \ \ ({-t \over 2} -3), \ \ \ \ (6t^2 + 24), 
 \ \ \  \ \ \ \ ({t^2 \over 3}+2t+20), \ ...$$
$$ where \ h_i = t*h_{i-2} + c_ih_{i-1} \\ and \ \ \ \ \ k_i = t*k_{i-2} + c_ik_{i-1} \\ i = 2,...,\infty$$
