I am using perturbation theory to solve a problem of the following form: $$ R(h,\theta)f(\theta) = h g(\theta) = 0 $$ where $h$ is small, and I assume $$ \theta = \sum_{i=0}^\infty \theta_i h^i $$ Further, I have expanded $R(h,\theta)$ as $$ R(h,\theta) = \sum_{i=0}^\infty R_i(\theta)h^i $$ and I expand $f$, $g$, and $R_i(\theta)$ as Taylor series \begin{align} f(\theta) &= f(\theta_0) + f'(\theta_0)(\theta-\theta_0) + \frac{1}{2}f''(\theta_0)(\theta-\theta_0)^2 + ...\\ g(\theta) &= g(\theta_0) + g'(\theta_0)(\theta-\theta_0) + \frac{1}{2}g''(\theta_0)(\theta-\theta_0)^2 + ...\\ R_i(\theta) &= R_i(\theta_0) + R_i'(\theta_0)(\theta-\theta_0) + \frac{1}{2}R_i''(\theta_0)(\theta-\theta_0)^2 + ... \end{align}
I would like Mathematica to list the terms of the equation I'm trying to solve at the various orders in $h$. For instance, at $O(1)$, we have just $$R_{00}f_0=0$$, where $R_{00}$ stands for $R_0(\theta_0)$ and $f_0$ stands for $f(\theta_0)$. Similarly, at $O(h)$ the term is $$(R_{00}'\theta_1 +R_{10})f_0 + R_{00}f_0'\theta_1+g_0 = 0$$ and at $O(h^2)$ we get $$ (R_{00}'\theta_1 + R_{10})f_0'\theta_1 + R_{00}(f_0'\theta_2 + \frac{1}{2}f_0''\theta_1^2) + g_0'\theta_1 = 0 $$ These terms are obviously becoming more complicated at the higher orders, and the chances of me making a mistake in finding these terms increases probably exponentially with the order of $h$. So I'm wondering, how can I get Mathematica to display these terms for me? I have tried already: here is what I'm using:
n = 5; θh := θ0 + Sum[θ[i] h^i, {i, 1, n}] + O[h]^(n + 1); rSymb[θ_] = Sum[rrExp[i, θ] h^i, {i, 0, n}] + O[h]^(n + 1); rrExp[i_, θ_] := Series[rr[i][θ], {h, 0, n}]; fSymb[θ_] := Series[ff[θ], {h, 0, n}]; gSymb[θ_] := Series[gg[θ], {h, 0, n}]; ySymb[θ_] := rSymb[θ] fSymb[θ] + h*gSymb[θ] rrExp is the Taylor expansion of the individual $R_i(\theta)$, fSymb and gSymb are the Taylor series expansions of $f$ and $g$, and ySymb is supposed to be the series consisting of the terms I want Mathematica to compute/display for me. However, there is a problem with rSymb. I don't know how to get it to properly regroup the rrExp terms by powers of $h$. I'm looking for rSymb to give me a series of the form \begin{multline} R(h,\theta) = R_{00} + h(R_{00}'\theta_1 + R_{10})\\ + h^2(R_{00}'\theta_2 + \frac{1}{2}R_{00}''\theta_1^2 + R_{10}'\theta_1 + R_{20})\\ +h^3(R_{00}'\theta_3 + R_{00}''\theta_1\theta_2 +\frac{1}{6} R_{00}'''\theta_1^3+R_{10}'\theta_2 + \frac{1}{2}R_{10}''\theta_1^2+R_{20}'\theta_1 + R_{30})\\ + O(h^4) \end{multline}
How can I modify my code to get Mathematica to properly find rSymb and then use that in ySymb to automatically determine the terms at arbitrary order in h?
Collectto the result ofrSymb? $\endgroup$