This article gives an expression for $\ln D$:
$$\ln D=\left( \begin{array}{cccccc} \psi (1) & 0 & 0 & 0 & 0 & . \\ 0 & \psi (2) & 0 & 0 & 0 & . \\ 0 & 0 & \psi (3) & 0 & 0 & . \\ 0 & 0 & 0 & \psi (4) & 0 & . \\ 0 & 0 & 0 & 0 & \psi (5) & . \\ . & . & . & . & . & . \\ \end{array} \right)-\ln x$$
I derived the same result using the following Mathematica code, starting from the first principles:
Table[DSolveValue[{f'[x] + s f[x] == x^n, f[0] == 0}, f[x], x] // FullSimplify, {n, 0, 5}]
Integrate[%, s] // FullSimplify
Limit[%, s -> 0, Direction -> "FromAbove"] // FullSimplify
Table[Coefficient[%, x, k], {k, 0, 5}] // Expand // MatrixForm
But when I try to exponentiate it, I do not obtain $D$. What's wrong?
