# Finding ODE series solution coefficients

I am trying to solve an ODE by subbing in a series form and then looking individually at the coefficients of different powers of the variable. I'm looking at a general form of equation:

$$\frac{\mathrm{d}^2w}{\mathrm{d}z^2} + f_1(z) \frac{\mathrm{d}w}{\mathrm{d}z} + f_0(z)=0$$

I want to fill in the solution form $$w(z) = e^{\lambda z} z^\mu \sum_{s=0}^\infty a_s z^{-s}$$ and then collect coefficients of the various powers of $$z$$. However when I try to do this the series terms are not being multiplied with any $$z$$ terms which are multiplying the series from the outside. A simpler example of the problem I'm having is when I input say:

w[z_] = E^(λ z)z^μ  Sum[a[s] z^-s, {s, 0, Infinity}]
Coefficient[z^5 w[z], z^μ]

The output of this is to give the coefficient as: $$e^{\lambda z} z^5\sum_{s=0}^\infty a_s z^{-s}$$

What I want is for the $$z^5$$ to multiply the $$z^\mu$$ term and also term-by-term with the series terms, and then give me just the coefficient of the resulting $$z^\mu$$ term. So the only corresponding term in the sum is the one of the form $$e^{\lambda z}z^5z^\mu a_5 z^{-5}$$, and therefore the coefficient I should be getting is just: $$e^{\lambda z} a_{5}$$

Can I not do this type of manipulation or have I just done something wrong here?

• I've tried to get mathematica to handle frobenius method solutions to ODEs in the straightforward, "code it like you would do it by hand" sort of way and I haven't been able to make it work to any satisfactory extent. My failure aside, you might want to see if the symbol AsymptoticDSolveValue is sufficient for your use-case. Depending on the specific form of your $f_1$, $f_0$ it may or may not suffice. Dec 23, 2019 at 23:38

Let me show you an example of series solution of ODE with given f0 and f1 for your differential equation.

{f0[z_] = Sin[z], f1[z_] = 1};

deq = Derivative[2][w][z] + f1[z]*Derivative[1][w][z] + f0[z] == 0

For comparison the solution with DSolve at initial conditions.

wdsol = w /. First@DSolve[deq && w[0] == 4 && w'[0] == 3, w, z]

(*   Function[{z}, 1/2 E^-z (-5 + 12 E^z + E^z Cos[z] + E^z Sin[z])]   *)

Generate a series expression. Don't omitt O[z]^smax+1. LogicalExpand generates equations for each power of z. Solve to get the coefficients a[s]. Get the last coefficients with initial conditions.

w[z_] = Sum[a[s] z^s, {s, 0, 15}] + O[z]^16

le = LogicalExpand[deq]

sol1 = Solve[le, Table[a[s], {s, 2, 15}]]

ww[z_] = Normal[w[z] /. First@sol1]

sol2 = Solve[{ww[0] == 4, ww'[0] == 3}, {a[0], a[1]}]

wle[z_] = ww[z] /. First@sol2

(*   4 + 3 z - (3 z^2)/2 + z^3/3 - z^4/12 + z^5/40 - z^6/240 + z^7/2520 - z^8/20160 \
+ z^9/120960 - z^10/1209600 + z^11/19958400 - z^12/239500800 + \
z^13/2075673600 - z^14/29059430400 + z^15/653837184000   *)

Plot[{wle[z], wdsol[z]}, {z, 0, 7}, PlotStyle -> {Blue, Red}]