Is there something faster than AsymptoticDSolveValue for getting the terms in a power series solution?

Question

Here's an example of a differential equation which Mathematica 13.1 just returns without solving

DSolve[{y''[x] + (5 x^3 + 3 x^2 + 5 x + 1) y'[x] + x y[x] == 0, y[0] == 1, y'[0] == 2}, y[x], x]

This specific differential equation isn't so important, it's just an example. If you want a power series solution around $x = 0$ you can do something like

AsymptoticDSolveValue[{y''[x] + (5 x^3 + 3 x^2 + 5 x + 1) y'[x] + x y[x] == 0, y[0] == 1, y'[0] == 2}, y[x], {x, 0, 1000}];

The problem is that I don't need 1,000 terms, I need more like 1,000,000 and that's going to take ages. I assume - perhaps wrongly - that if you painstakingly derived the recurrence relation for the coefficients in the power series solution, a command like RecurrenceTable[] would get the terms much faster, especially if we only need floats and not the exact values. But is there an automatic way to do this?

Answer 1

You can't obtain the recurrence formula with V 13.1

Looking at "Q&A with Calculus Developers: Live with the R&D team" Nov 30, 2022 at https://www.twitch.tv/videos/1666674721 for V 13.2, I did not see this feature being there (see around 7 minutes frame). Not even Maple can do this for series solutions. i.e. give explicit recurrence formula used for series solution.

I have a function myself that does this (but still in development) and it only handle linear 1st and second order odes'. For your ode, this is the recurrence formula it found

$n=0$ gives $$ a_{2} = -\frac{a_{1}}{2} $$ $n=1$ gives $$ a_{3} = -\frac{a_{0}}{6}-\frac{2 a_{1}}{3} $$ $n=2$ gives $$ a_{4} = \frac{a_{0}}{24}+\frac{a_{1}}{4} $$ For $3\le n$, the recurrence equation is $$ a_{n +2}= -\frac{5 n a_{n}}{\left(n +2\right) \left(n +1\right)}-\frac{\left(5 n -10\right) a_{n -2}}{\left(n +2\right) \left(n +1\right)}-\frac{\left(3 n -2\right) a_{n -1}}{\left(n +2\right) \left(n +1\right)}-\frac{a_{n +1}}{n +2} $$

Now you are able to generate very fast all terms you want. Here is the Mathematica code

a[2] = -a[1]/2
a[3] = -a[0]/6 - 2 a[1]/3
a[4] = a[0]/24 + a[1]/4
Table[a[n + 2] = -5*n*a[n]/((n + 2)*(n + 1)) - (5*n - 10)*
     a[n - 2]/((n + 2)*(n + 1)) - (3*n - 2)*
     a[n - 1]/((n + 2)*(n + 1)) - a[n + 1]/(n + 2), {n, 3, 10}];
sol = Sum[a[n]*x^n, {n, 0, 10}];
Collect[sol, {a[0], a[1]}]

Compare to

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

Note that I generates all the $a_n$ using normal Table command. You could try to use RecurrenceTable instead to see if performance is better.

You can now solve for $a_0,a_1$ easily given the initial conditions.

Answer 2

You can obtain the coefficients using SeriesCoefficient on the DifferentialRoot of the differential equation

f = DifferentialRoot[
  Function[{y, 
    x}, {y''[x] + (5 x^3 + 3 x^2 + 5 x + 1) y'[x] + x y[x] == 0, 
    y[0] == 1, y'[0] == 2}]]

rec = SeriesCoefficient[f[x], {x, 0, n}][[1, 1, 1, 0, 1]][y, n]

(* {5 n y[n] + (4 + 3 n) y[1 + n] + (10 + 5 n) y[2 + n] + (3 + n) y[ 3 + n] + (3 + n) (4 + n) y[4 + n] == 0, y[0] == 1, y[1] == 2, y[2] == -1, y[3] == -(3/2)} *)

In LaTeX :

$$ \left\{5 n y(n)+(3 n+4) y(n+1)+(5 n+10) y(n+2)+(n+3) y(n+3)+(n+3) (n+4) y(n+4)=0,y(0)=1,y(1)=2,y(2)=-1,y(3)=-\frac{3}{2}\right\} $$

To compare with @Nasser 's recurrence we can use:

Solve[rec[[1]] /. n -> n - 2, y[2 + n]][[1, 1]] // Simplify

(* y[2 + n] -> -( 1/((1 + n) (2 + n)))(5 (-2 + n) y[-2 + n] + (-2 + 3 n) y[-1 + n] + 5 n y[n] + y[1 + n] + n y[1 + n]) *)

In LaTeX

$$ y(n+2)\to -\frac{5 (n-2) y(n-2)+(3 n-2) y(n-1)+5 n y(n)+n y(n+1)+y(n+1)}{(n+1) (n+2)} $$