# Solving an ODE in power series

How do I find a series solution to an ODE? I do not mean taking the Taylor series of an exact solution; I want to solve nasty nonlinear differential equations locally via plug and chug. Surely, that is a built-in option in Mathematica.

• Yes, there's a way. Do you have an example problem? Commented May 17, 2013 at 3:45
• Somewhat surprisingly, the procedure you want isn't built-in, but one could certainly write a routine for this task. Commented May 17, 2013 at 3:46
• @J.M. Not true! It's buried in the Holonomic context Commented May 17, 2013 at 3:50
• @Mark, I see. I had something like this in mind... still, if memory serves, the functions there can only do linear ODEs, not things like $y^\prime=1+y^2$. Commented May 17, 2013 at 3:55
• @Mark, after playing with some of those functions for a while, I'm not sure if using them is any better than forming a DifferentialRoot[] and then using Series[] for expansion (altho I am sure those functions are used under the hood with what I'm proposing). Commented May 17, 2013 at 5:06

This method of developing a truncated solution can be done as below. I illustrate with an example that DSolve does not seem much to like.

ode = x''[t] - t*x'[t] + Sin[t] == 0;
initconds = {x'[0] == 1, x[0] == 0};


We create a differentail operator to create this ode.

odeOperator = D[#, {t, 2}] - t*D[#, t] + Sin[t] &;


Now set up our Taylor series as a symbolic expansion using derivatives of x evaluated at the origin. I use an order of 15 but that is something one would probably make as an argument to a function, if automating all this.

xx = Series[x[t], {t, 0, 15}];


Next apply the differential operator and add the initial conditions. Then find a solution that makes all powers of t vanish.

soln = SolveAlways[Join[{odeOperator[xx] == 0}, initconds], t];


Let's look at the Taylor polynomial.

truncatedSol = Normal[xx /. soln[[1]]]

(* Out[500]= t + t^5/120 + t^7/1260 + (29 t^9)/362880 + (
13 t^11)/1995840 + (2861 t^13)/6227020800 + (4649 t^15)/163459296000 *)


To assess how accurate it might be I will compare with NDSolve.

approxSol = NDSolve[Join[{ode}, initconds], x[t], {t, 0, 4}];

Plot[{truncatedSol, x[t] /. approxSol[[1]]}, {t, 0, 4}]


• +1 This is a good start. Power series solutions, though, are frequently used to obtain recursion equations for the coefficients (of any solution that might be analytic within a neighborhood of the point of expansion). It would be nice, then, to have a function that outputs these equations (given a differential operator as input), rather than just obtaining an approximate solution with a limited radius of accuracy. In order to analyze singular points, it would also be useful to consider slightly more general series of the form $z^\alpha(a_0+a_1z+a_2z^2+\cdots)$ for non-integral $\alpha$. Commented May 17, 2013 at 18:07
• Thank you for your answer; it is very clear. My equation is very similar to the one you use; its nonlinear term is of the form sin(2f). I am more familiar with Maple which has the basic command dsolve(equation, series) and I was sure that MMA must have such a command also.
– rick
Commented May 18, 2013 at 2:02
• @rick, that's a bit unclear; is the thing inside the sine the dependent or independent variable? Commented May 18, 2013 at 15:03
• My equation has the term sin (2y(x)) so it is far from linear. When looking at the documentation for NDSolve I did not see see where one might specify the integration method to use-I was searching for a Taylor series option.
– rick
Commented May 19, 2013 at 15:25
• This's a great answer, but I don't know why we can't compute the following differential equation case:ode = x[t]*x''[t] - 1 - x'[t]^2 == 0; initconds = {x[0] == 1, x'[0] == 0};. Commented Aug 23, 2020 at 3:25

J.M. in a comment above alluded to this possibility. Using Series and DifferentialRoot seems like a great way to expand linear differential equations in series. It seems worth putting down in an answer due to its simplicity:

lde = {y''[x] - x*y'[x] + Sin[x] == 0, y'[0] == 1, y[0] == 0};
Series[DifferentialRoot[Function @@ {{y, x}, lde}][x], {x, 0, 15}]

(* x + x^5/120 + x^7/1260 + (29 x^9)/362880 + (13 x^11)/1995840 +
(2861 x^13)/6227020800 + (4649 x^15)/163459296000 + O[x]^16 *)


Here's a way to develop a series solution for a quasilinear ODE with analytic coefficients, $y^{(n)}=f(x,y,y',\dots,y^{(n-1)})$. The main idea to fold Solve over the series expansion of the ODE solve for the derivatives of $y$, something like this:

Fold[
Flatten[{#1, Solve[#2 == 0 /. #1]}] &,
ics,                                               (* initial conditions *)
Series[ode /. Equal -> Subtract, {x, x0, n}][[3]]  (* series coefficients *)
]


Getting this idea to work for others without their having tweaking it takes some work. I've handled most of the things I've run into in the code below, but there are probably some things I have not tested. A malformed ODE may fail, probably ungracefully. Suggestions and comments are welcome if you run into other problems.

Clear[seriesDSolve];
seriesDSolve[ode_, y_, iter : {x_, x0_, n_}, ics_: {}] :=
Module[{dorder, coeff},
dorder = Max[0, Cases[ode, Derivative[m_][y][x] :> m, Infinity]];
(coeff = Fold[
Flatten[{#1,
Solve[#2 == 0 /. #1,
Union@Cases[#2 /. #1, Derivative[m_][y][x0] /; m >= dorder,
Infinity]]}] &,
ics,
Table[
SeriesCoefficient[ode /. Equal -> Subtract, {x, x0, k}],
{k, 0, Max[n - dorder, 0]}]
];
Series[y[x], iter] /. coeff) /; dorder > 0
]


Example (constructed from a comment by the OP):

seriesDSolve[
y''[x] + y'[x] == Sin[2 y[x]],
y, {x, 0, 3},
{y[0] -> a, y'[0] -> b}         (* symbolic initial values *)
]


(Looking at it, I suppose it would be nice if the initial condition could be specified just like in DSolve. That would take some parsing and this approach evolved according to in-house needs -- that is to say, mine.)

• As we know, any Nth order ordinary differential equation can be recast as a coupled set of N ﬁrst order ODEs, but not every set of N ﬁrst order ODEs can necessarily be written as a neat Nth order ODE. So , if there are a system of coupled ODEs, how can we use seriesDSolve to obtain the series solution? Commented Nov 9, 2016 at 15:41
• @user10709 You'd have to change the code to handle a vector of dependent variables instead of just a single (scalar) y. That includes handling their different derivatives by adapting the pattern Derivative[m_][y][x0]. The Solve command would have to solve a system for the next higher-order derivatives of the variables at each stage. That seems feasible, but easier if it's a first order system of the form $dx/dt = F(x, t)$ (for vector $x$). As you say, a system could be put in such form. Commented Nov 9, 2016 at 18:55
• @MichaelE2, This is very useful code. I used it to determine series solution for second order nonlinear differential equation. Can you extend it to four coupled second order differential equation? Commented Jan 30, 2021 at 0:45

On Mathematica 11.3 using a new function AsymptoticDSolveValue

Linear ODE:

AsymptoticDSolveValue[{y''[x] - x*y'[x] + Sin[x] == 0, y'[0] == 1,
y[0] == 0}, y[x], {x, 0, 10}]

(* x + x^5/120 + x^7/1260 + (29 x^9)/362880  *)


Quasilinear ODE:

AsymptoticDSolveValue[{y''[x] + y'[x] == Sin[2*y[x]], y[0] == a,
y'[0] == b}, y[x], {x, 0, 5}]

(* a + b x + 1/6 x^3 (b + 2 b Cos[2 a] - Sin[2 a]) +
1/2 x^2 (-b + Sin[2 a]) +
1/24 x^4 (-b - 4 b Cos[2 a] + Sin[2 a] - 4 b^2 Sin[2 a] +
2 Cos[2 a] Sin[2 a]) +
1/120 x^5 (b + 6 b Cos[2 a] - 8 b^3 Cos[2 a] + 4 b Cos[2 a]^2 -
Sin[2 a] + 16 b^2 Sin[2 a] - 4 Cos[2 a] Sin[2 a] -
12 b Sin[2 a]^2)*)


Nonlinear ODE:

AsymptoticDSolveValue[{y'[x] == 1 + y[x]^2, y[0] == 0}, y[x], {x, 0, 10}]

(* x + x^3/3 + (2 x^5)/15 + (17 x^7)/315 + (62 x^9)/2835 *)
`