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.
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] 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.
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]
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$, and checked it against your paper. 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$.
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/2+1): #there was a typo here
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;
[Edited to fix the above code.]
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.
[Edit]
The generic coefficients in the recursive formula only depend on the polynomials $P_i$. I'll explain when each one is homogeneous of degree $d_i$. The inhomogeneous case is not nearly as clean, so I'd introduce a new variable with substitution polynomial $1$ to reduce it to this case. To find the coefficient of $f_{n_1,\ldots,n_k}$, write $T=\sum_i n_i$ for the total degree. Then in the recursive formula, the smaller coefficient $f_{m_1,\ldots,m_k}$ can only appear when $T = \sum d_im_i$. To find its multinomial and constant parts, write equations for each subscript in terms of the polynomials $P_i$. If $P_i$ has $j_i$ many terms, you will need $j_i$ indices, say $\tau_{i,1}, \ldots, \tau_{i,j_i}$ which sum to $m_i$. These correspond to expanding $P_i^{m_i}$ via the multinomial theorem, which will give you $\binom{m_i}{\tau_{i,1}, \ldots, \tau_{i,j_i}}$ multiplied by terms $c_{i,\alpha_1, \ldots, \alpha_k}^{\tau_{i,h}}$, where the $c_i$s are the coefficients of and the $\alpha_1, \ldots, \alpha_k$ are the exponents of term $h$ in $P_i$. Yikes!
So now we know what the coefficients in the recursive formula should look like, but we need to sum only over the indices $\tau$ which send us to $f_{n_1,\ldots,n_k}$. Each variable gives us one equation, like $n_1 = \sum_{i,h} \tau_{i, h}\alpha_{1,h}$ where the $\alpha$ terms are the exponents again, and $h$ runs from $1$ to $j_i$ for each $i$. These, along with the defining equations $\sum_h\tau_{i,h}=m_i$ and $T = \sum d_im_i$ give us $2k+1$ linear equations ($1$ for each $n$, one for each $m$, and one for the total degree) with $k+1$ dependencies, so you can reduce it to just $k$ equations, and the sum is taken over all non-negative integer solutions to those.
So, for example, the coefficients of your function $G$ are
$$ g_{n_1,n_2,n_3} = \sum \binom{m_1}{\tau_{1,1}, \tau_{1,2}, \tau_{1,3}}\binom{m_3}{\tau_{3,1}, \tau_{3,2}}2^{\tau_{1,2}}g_{m_1,m_2,m_3} $$
Where the sum runs over all non-negative $\tau, m$ such that:
$$ 2m_1 + m_2 + 2m_3 = n_1 + n_2 + n_3 $$ $$ m_1 = \tau_{1,1} + \tau_{1,2} + \tau_{1,3} $$ $$ m_2 = \tau_{2,2} $$ $$ m_3 = \tau_{3,1} + \tau_{3,2} $$ $$ 2\tau_{1,1} + \tau_{1,2} + \tau_{2,1} + 2\tau_{3,1} + \tau_{3,2} = n_1 $$ $$ \tau_{1,2} + \tau_{1,3} + \tau_{3,2} = n_2 \\ $$ $$ \tau_{1,3} = n_3 $$ Sorry I couldn't get the align environment to work, but you get the picture. You can reduce it to just three equations, one of which is constant, which looks a lot nicer, but that's the generic way to do it.