You can get recurrence relations for the coefficients by considering, for each monomial, which lower degree monomials will contribute to its coefficient. For example, suppose we want to find the coefficients of $S$ without iteration. Iteration here is particularly bad: the degree of $S_{k+1}$ will be $3^k$, so it will contain monomials of degree $3^{k-1}+1$ which did not appear in $S_k$. In the limit, these will contribute to monomials of degree at least $2(3^{k-1}+1)$, so only about two thirds of the coefficients in $S_{k+1}$ will be accurate. Instead, consider which smaller terms will contribute to the coefficient of $x^n$. For any $m<n$, $x^m$ will contribute if $x^n$ appears in $$ (x^2 + x^3)^m = \sum_{i=0}^m \binom{m}{i}x^{3i}x^{2m-2i} $$ Thus, we want $n=2m+i$ for some $0 \leq i \leq m$. Write $j=m-i$, and $c_k$ for the coefficient of $x^k$ in $S$. Then to find $c_n$, we will sum $c_m \binom{m}{i}$ over all pairs of non-negative integers $(i,j)$ such that $2i+3j=n$ and $i+j=m$, i.e. $$ c_n = \sum_{2i+3j=n}\binom{i+j}{i}c_{i+j} $$ So to get $c_n$, we only need the coefficients with indices in the range $[n/3, n/2]$. Here is a bit of sage code which calculates the coefficients up to $k$, then tells us the last one. It finds the first thousand terms listed at OEIS A014535 in about 20 seconds for me. <pre><code>k=1000; #number of coefficients you want S=[0,1]; #list of coefficients, initialized with the first two values for n in range (2, k+1): count=0; #keeps track of the partial sums of the coefficient for i in range (0,n/2+1): for j in range((n-2*i)/3,n/3+1): if (2*i +3*j)==n: count=count+binomial(i+j,i)*S[i+j] S.append(count) S[k] </code></pre> Applying this to the problems with more variables, we find the coefficient $b_{m,n}$ of $x^my^n$ in $B$ is given by $$ b_{m,n} = 2^n \sum_{2i+j = n+m}\binom{i}{n}b_{i,j} $$ In particular, this tells us something we already know from finding $B$ by iteration: the coefficient is always $0$ when $n>m$, since in that case, $2i \leq 2i+j < 2n$, so all of the binomial coefficients in the sum are $0$. This lets us take a little shortcut in the code and skip all of those terms. <pre><code> import numpy k=500; B=numpy.zeros((k+1,k+1), dtype=object); #an array to hold the coefficients, with object type to give arbitrary precision integers B[1,0]=1; #initial values from the polynomial B[1,1]=2; for m in range(2,k+1): for n in range(0,m+1): count=0; for i in range (0, (m+n)/2 + 1): for j in range (0, i+1): if 2*i+j==m+n: count=count+(2^n)*(binomial(i, n))*B[i,j] B[m,n]=count; B[k,0] </code></pre> This gives us the last term listed in the OEIS A006265 (as well as all of the previous terms, and all of the coefficients with multiples of $y$), though it takes quite a bit longer (for me, less than 15 minutes but more than 3). Finally, I wrote up some code for your function $F$, but I don't have any place to check the values, so although the reasoning about the recurrence is the same, it may have some errors. The coefficients $f_{k,m,n}$ of $x^ky^mz^n$ are given by $$ f_{k,m,n} = \sum_{(a,b,c)}\binom{a}{n}\binom{a-n}{m-c-n}f_{a,b,c}; $$ where the sum is over all non-negative integers with $2a+b+2c=k+m+n$, so each one is determined by coefficients with strictly smaller total degree (since the monomials $y^m$ never show up). I split the multinomial into two binomials to make some of the restrictions more obvious: for example, if the degree of $z$ is more than the degree of $y$, then the coefficient is 0, because the second binomial always vanishes. This is clear from the given polynomials: the only monomial with a $z$ also contains $y$. Again, I'm not sure if I implemented this properly, but I'll give some results near $50$ and you can tell me if they look right: <pre><code> j=60; F=numpy.zeros((j,j,j), dtype=object); F[1,0,0]=1; F[2,0,0]=1; F[1,1,0]=1; F[0,1,1]=1; for t in range(3,j): #t is the total degree of the monomial for k in range (0, t+1): #k is the degree of x for m in range (0, t+1-k): #m is the degree of y n=t-k-m; #n is the degree of z if m >=n: count=0; for c in range(0, m-n+1): for a in range (max(n,m-n-c), (t+1)/2): b=t-2*a - 2*c; if b >=0: count=count+binomial(a,n)*binomial(a-n,m-c-n)*F[a,b,c]; F[k,m,n]=count; for r in range(54, j): F[r,0,0] </code></pre> So from 54-59 inclusive, I got: <code> 60346891 49793458 38341939 31272992 42718324 112825552 </code> I picked that string because they follow some coefficients which about the same size as the one on $x^{54}$, then decrease, then increase again. Seems unlikely to accidentally agree with your function, and may indicate a similar oscillation. A few other remarks: You will be able to find recurrence relations depending on coefficients of lower total degree for your other functions, as well as any where all polynomials have minimal degree $2$, or, those where if $P_i$ has a term of degree $1$, then $x_i^n$ never appears. This might tell you what you need to know about d-finiteness, though I don't know all the intricacies in that topic. Mill, Pippenger, Rosenberg and Snyder [SICOMP, 1979] gave the recurrence relation for $S$ above, arguing combinatorially, and gave asymptotics for it in a manner similar to the central limit theorem, by estimating the binomial coefficients. Such an argument may also give you asymptotics in the other cases, though you may need to pass to higher dimensional central limit arguments to deal with the multinomials.