I'm trying to solve an equation of the form $$ R(\theta)f(\theta)+hg(\theta) = 0 $$ for small $h$, where $R$, $f$, and $g$ are functions. I've assumed a power series expansion for $\theta$ in terms of $h$: $$ \theta = \sum_{i=0}^\infty \theta_i h^i $$ and I'm wondering how I can solve this in Mathematica in sort of an automated way. One complication is that I don't have an explicit expression for $R(\theta)$. Instead, I have $R'(\theta) = m(\theta)R(\theta)$ for some function $m$, and I know $R(\theta_0)\equiv R_0 \neq 0$. I do, however, have explicit expressions for $f(\theta)$ and $g(\theta)$. Assuming the following expansions $$ R(\theta) = R_0 +R'(\theta_0)(\theta-\theta_0) + \frac{1}{2}R''(\theta_0)(\theta-\theta_0)^2 + ...\\ f(\theta) = f_0 + f'(\theta_0)(\theta-\theta_0) + \frac{1}{2}f''(\theta_0)(\theta-\theta_0)^2 + ...\\ g(\theta) = g_0 + g'(\theta_0)(\theta-\theta_0) + \frac{1}{2}g''(\theta_0)(\theta-\theta_0)^2 + ... $$ where $f_0 \equiv f(\theta_0)$ and $g_0 \equiv g(\theta_0)$. I can find the first few terms by hand. For instance, at $O(1)$: $$ R_0f_0 = 0 \implies f_0=0 $$ This gives me an expression for $\theta_0$. Then, at $O(h)$, I obtain $$ (R'(\theta_0)f_0 + R_0f'(\theta_0))(\theta-\theta_0) + g_0 = 0 $$ Since we already determined $f_0 = 0$, this becomes $$ R_0f'(\theta_0)(\theta-\theta_0) + g_0 = 0 $$ Then, using the expansion for $\theta$, and remembering this is the $O(h)$ term, we obtain $$R_0f'(\theta_0)\theta_1 + g_0 = 0\implies \theta_1 = -\frac{g_0}{R_0 f'(\theta_0)}$$ and I can continue this. The question is how can I do this in Mathematica. I have tried several versions of the following:
θw[n_] := Sum[θ[i] h^(i - 1), {i, 1, n}] + O[h]^n;
f[θ_, n_] := Series[f[θw[n]], {h, θ0, n}]
g[θ_, n_] := Series[g[θw[n]], {h, θ0, n}]
R[θ_, n_] := Series[R[θw[n]], {h, θ0, n}]
j[θ_, n_] := R[θ, n] f[θ, n] + h*g[θ, n]
Solve[j[θ, 2] == 0]
{{g[θ[1]] -> -R[θ[1]] θ[2]
Derivative[1][f][θ[1]],
f[θ[1]] -> 0}, {g[θ[1]] -> -f[θ[1]] θ[2] Derivative[1][
R][θ[1]], R[θ[1]] -> 0}}
The part in the curly braces is the output of the Solve command. The $n$ in the functions is to set the order in $h$ to which the expansions are taken. This seems close to working. It does give me something like the expressions I derived. However, I would like for it to solve for the $\theta_i$s automatically, and I would like for it to ignore the $R_0=0$ solution, because I'm not interested in that one. A further complication is that for higher order corrections, I will need to implicitly differentiate my expression for $R$: $$R'(\theta) = m(\theta)R(\theta)$$ Further, is it possible to automatically use the information from lower orders to find explicit expressions for the higher order $\theta_i$s?