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?

  • $\begingroup$ 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. $\endgroup$
    – Diffycue
    Commented Dec 23, 2019 at 23:38

1 Answer 1


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}]

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.