Although I don't have a "final" answer, this suggestion may help.

For polynomials or for analytic functions (having a power series representation with nonzero radius of convergence) I'd employ the concept of Carleman-matrices.

Assume a vectorfunction $V(x) = [1,x,x^2,x^3,...]$ as rowvector and $F$ as carleman-matrix (transposed) for your function $f(x)$ and *I* for the identity-matrix then we could in principle write

$\qquad V(a_1) \cdot I = V(a_1)$

$\qquad V(a_1) \cdot F = V(f(a_1)) $

but in

$\qquad V(a_1) \cdot (I+F) = V(a)+ V(f(a_1)) \ne V(a_2)$

the sum of two *V()*-vectors is not a *V()*-vector.

Instead we define first the Carleman-matrix $G$ for the function $g(x)=x+f(x)$

Then we can iterate:

$\qquad V(a_0) = V(a_0) \cdot I\\
\qquad V(a_1) = V(a_0) \cdot G \\
\qquad V(a_2) = V(a_1) \cdot G = V(a_0) \cdot G^2 \\
\qquad \cdots \\
\qquad V(a_k) = V(a_0) \cdot G^k \\
$

as long as taking the *k*'th power $G^k $ makes sense (requires only convergent (or as generalization for certain divergent cases for instance Euler-summable) series).

If $G$ is **triangular**, the formal power series for your iterated expression $a_k$ can exactly be given to any power (even for fractional powers!) and with your initial value $a_o$ it might be convergent up to some highest power of *G* . For instance, the transposed Carleman-matrix for the function $ f(x) = \ln(1+x)$ *is* triangular, and also that for $g(x) = x + \ln(1+x) $, however the range of convergence of the formal power series decreases with the iterations...

If $G$ is triangular and $g(x)' = 1$ we can express the power series for the k'th iteration with coefficients, which are polynomials in *k* and are thus especially easy to compute.

If $G$ is triangular and $g(x)' \notin (0,1) $ we can apply matrix-diagonalization, which gives again exact coefficents for the power series, but are likely more complicated.

*(See two examples below)*

For **non-triangular** Carleman-matrices (where $g(0) \ne 0$) this is even more delicate and exceeds the focus of this answer...

Two examples.
For

$f(x) = \ln(1+x)$ the power-series for the *h*-fold iterate (

$f°^h(x)$) begins with

$ \small \begin{eqnarray} f°^h(x) =&
x \cdot &(1) \\
&+x^2 \cdot &-( 1/2 h) \\
&+x^3 \cdot & (1/12 h&+1/4 h^2 ) \\
&+x^4 \cdot & -( 1/48 h &+ 5/48 h^2& +1/8 h^3) \\
&+x^5 \cdot & ( 1/180 h& + 1/24 h^2& + 13/144 h^3&+1/16 h^4 ) \\
&+ O(x^6) \end{eqnarray}
$

and for $g(x)= x + \ln(1+x) $ the powerseries for the *h*-fold iterate begins with

$ \small \begin{eqnarray} g°^h(x) =&
x \cdot &(1 u) \\
&+x^2 \cdot &(1/4 u-1/4 u^2) \\
& +x^3 \cdot &(1/36 u-1/8 u^2+7/72 u^3) \\
& +x^4 \cdot &(1/672 u-17/576 u^2+7/96 u^3-181/4032 u^4) \\
& +x^5 \cdot &(11/75600 u-17/4032 u^2+91/3456 u^3-181/4032 u^4+13687/604800 u^5)\\
&+O(x^6) \\
\end{eqnarray} $

where we have to replace $u$ by $2^h$ and *h* is the iteration-height-parameter.