Remark: I give much background because it might significantly help to find an idea how to generate a solution
I'm analyzing properties of a certain infinite matrix $U$, for whose columns we have the generating functions $$f_c(x) = (\exp(x)-1)^c \qquad \qquad c=-\infty...+\infty $$, assuming $c$ from $-\infty$ to $+\infty$ (This is the type of Carleman-matrices extended to two-way-infinite size which Eri Jabotinsky has discussed in one article I've found in some archive...)
Clearly the Stirling numbers 2nd kind are involved in the part of nonnegative column-indexes, but in the columns with negative index we have Zeta-values (or: Bernoulli-numbers) due to the generating functions having $\exp(x)-1$ to the negative power.
I succeded with finding that the generating functions for the rows are $$ g_r(x) = {1 \over 1+x}{1 \over \log(1+1/x)^{1+r}} \qquad \qquad c=-\infty...+\infty$$ for rowindex r from $-\infty$ to $+\infty$ .
(We can find the coefficients easier when we write $$ g^*_r(t)={t \over 1+t}{1 \over \log(1+t)^{1+r}} $$ get the Laurent-series for this and replace the powers of t by the negative powers of x .)
Now, because I'm looking at iterations of the exponential function I've the matrix $U$ to the second power, $U^2$ which has the following entries around its zero-indexes:
Because the 2nd power of the matrix gives the second iterate of the function in the columns, I know that the generating functions of the columns are now $ f_{2,c}(x) = (\exp(\exp(x)-1)-1)^c$
Q: My problem is now to find the generating functions of the rows. Simply writing the iterate of the rows-functions of the original version of $U$ seems to lead to nothing when I let Pari/GP do the expression. Does someone see how I could express them given the beginning of the sequence of coefficients?
For instance, at row $r=0$ whe have $$g_{2,0}(x) = 1 -1/x + 4/3/x^2-15/8/x^3+122/45/x^4 + \cdots = ??? $$ or
$$g^*_{2,0}(t) = 1 -1t + 4/3t^2-15/8t^3+122/45t^4 + \cdots = ??? $$
A meaningful suggestion for only one of the rows might already be helpful because the overall system should have a simple systematic derived from one key-idea as it is also for the non-iterated case.
Ahh, this looks very good...
By some trial & error I got the following scheme, which extends also nicely to higher powers of $U$ (or: iterates of $\exp(x)-1$). I got the following four families of generating functions for iteration-counts $h=0,1,2,3$ depending on rowindex $r$: $$ \small \begin{array} {} U^0:& g^*_{0,r}(t)&=& t \cdot{1 \over t^{1+r} } \qquad \qquad \text{note: $U^0$ is the identity-matrix}\\ U^1:& g^*_{1,r}(t)&=& t \cdot {1\over 1+t}\cdot{1 \over \log(1+t)^{1+r} } \\ U^2:& g^*_{2,r}(t)&=&t \cdot {1\over 1+t}{1 \over 1+\log(1+t)}\cdot{1 \over \log(1+\log(1+t))^{1+r}} \\ U^3:& g^*_{3,r}(t)&=& t \cdot {1\over 1+t}{1 \over 1+\log(1+t)}{1 \over 1+\log(1+\log(1+t))}\cdot {1 \over \log(1+\log(1+\log(1+t)))^{1+r}} \\ \end{array} $$ and I assume that this is the general pattern also for the higher families.
Denoting $\Lambda(t)=\log(1+t)$ and for the iterations $\Lambda^{°h}(t)$ where $\Lambda^{°0}(t)=t$ we seem to have $$ g^*_{h,r}(t) = \left(\prod_{k=0}^{h-1} {1 \over 1+\Lambda^{°k}(t)} \right) \cdot {t \over \Lambda^{°h}(t)^{1+r}} $$
Still open: an analytical justification/confirmation for this heuristic...
Update 10. Jun 2015: In this answer there is an analytical derivation valid for more general cases (functions of the type $f(x) = x + O(x^2)$ )