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
we find that the first 22 terms of nr(nr(nr(1+t))) are correct.
Your task is to write an efficient implementation of this method.