We can find a solution to our equation $C(t) =\Phi(C(t))$ more efficiently by using Newton-Raphson. Note that $\Phi$ is a power series (actually polynomial in $u$) and \[ \Phi(C_0+\de) \approx \Phi(C_0) +\Phi'(C_0)\de; \] this works provided $\de$ is divisible by a power of $t$. So if $C_0$ is an approximate solution to $\Phi(C(t))=C(t)$, then our aim is to choose $\de$ so that \[ C_0 + \de = \Phi(C_0) +\Phi'(C_0)\de, \] which implies that \[ \de = \frac{\Phi(C_0)-C_0}{1-\Phi'(C_0)}. \] Hence our proposed update formula is now \[ C_{n+1} = C_n + \frac{\Phi(C_n)-C_n}{1-\Phi'(C_n)} \] Here \[ \Phi'(C_n(t)) = 2tC_n(t). \]

We can run a quick test. With

phi = lambda u: 1+t*u^2
nr = lambda u: u + (phi(u)-u)/(1-2*t*u)

we find that the first 22 terms of nr(nr(nr(1+t))) are correct.

nr(nr(nr(1+t)))

1 + t + 2*t^2 + 5*t^3 + 14*t^4 + 42*t^5 + 132*t^6 + 
429*t^7 + 1430*t^8 + 4862*t^9 + 16796*t^10 + 58786*t^11 + 
208012*t^12 + 742900*t^13 + 2674440*t^14 + 9694845*t^15 + 
35357670*t^16 + 129644790*t^17 + 477638700*t^18 + 
1767263190*t^19 + 6564120420*t^20 + 24466267020*t^21 + O(t^22)

nr(nr(nr(1+t))).O(19) - Ct.O(19)

O(t^19)

Your task is to write an efficient implementation of this method.