# Why both DSolve and NDSolve are unable to solve a second-order differential equation?

I am trying to solve a recurrence relation using generating functions method: $$a_n=a_{n-1}+(n-1)a_{n-2}+(0.5n^2-1.5n+1)a_{n-3}$$

After some long calculations, I have arrived to this second-order differential equation: $$0.5 x^5 y''(x)+(2x^4+x^3)y'(x)+\left(x^3+x^2+x-1\right)y(x)+1=0$$

and these conditions: $$y(0)=1, y'(0)=1$$. $$y(x)$$ is the function that needs to be expanded as Taylor Series at $$x=0$$ to obtain the sequence from the coefficients. However, when I try to solve it both using DSolve and NDSolve, I have no luck. With DSolve it just returns the request itself:

$$\text{DSolve}\left[\left\{0.5 x^5 y''(x)+(2. x+1) x^3 y'(x)+\left(1. x^3+x^2+x-1\right)y(x)+1=0,y(0)=1,y'(0)=1\right\},y,x\right]$$

And with NDSolve I just receive errors and no equation:

Power::infy: Infinite expression 1/0.^5 encountered.
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0..


$$\text{NDSolve}\left[\left\{0.5 x^5 y''(x)+(2. x+1) x^3 y'(x)+\left(1. x^3+x^2+x-1\right)y(x)+1=0,y(0)=1,y'(0)=1\right\},y,\{x,0,1\}\right]$$

How could I resolve this problem?

Code:

Simplify[y[x] - (1 + x + 2 x^2)]
l = Expand[%]
Simplify[x (y[x] - (1 + x))]
r1 = Expand[%]
Simplify[x*D[x^2 (y[x] - 1), x] - x^2 (y[x] - 1)]
r2 = Expand[%]
Simplify[0.5 x*D[x*D[x^3*y[x], x], x] - 1.5 x*D[x^3*y[x], x] +
x^3*y[x]]
r3 = Expand[%]
eq = FullSimplify[r1 + r2 + r3 - l]
DSolve[{eq == 0, y[0] == 1}, y, x]
NDSolve[{eq == 0, y[0] == 1, y'[0] == 1}, y, {x, 0, 1}]

• Please post code that can be readily copy-pasted by anyone who wants to try and answer your question. Oct 14, 2020 at 17:44
• NDSolve fails here, because the ODE is singular at x = 0. Try integrating from a slightly larger value of x, but even then expect some difficulties. Oct 14, 2020 at 18:47
• RecurrenceTable[{a[n] == a[n - 1] + (n - 1) a[n - 2] + (n^2 - 3 n + 2)/2 a[n - 3], a[1] == 1, a[2] == 1, a[3] == 2}, a, {n, 11}] yields {1, 1, 2, 8, 22, 82, 334, 1370, 6338, 30692, 155722}. Oct 15, 2020 at 11:30
• I'm starting from 0, a[0]==1, a[1]==1, a[2]==2 Oct 16, 2020 at 11:33

Error Corrected.

Because neither DSolve nor AsymptoticDSolveValue produced a solution, I expanded the ODE itself and computed the coefficients of each term.

nmax = 10;
Normal@Series[x^5  y''[x]/2 + (2 x + 1) x^3 y'[x] + (x^3 + x^2 + x - 1) y[x] + 1,
{x, 0, nmax}]
CoefficientList[%, x] /. Derivative[n_][y][0] -> c[n] /. y[0] -> c[0]
Flatten@Solve[Thread[% == 0], Array[c, nmax + 1, 0]] // Values
%/(Range[0, nmax]!)

(* {1, 1, 2, 5, 14, 46, 166, 652, 2780, 12644, 61136} *)


When I rationalize eq, AsymptoticDSolveValue returns a solution for me:

eq = FullSimplify[r1 + r2 + r3 - l]

1 + (-1 + x + x^2 + x^3) y[x] + x^3 (1 + 2 x) y'[x] + 1/2 x^5 y''[x]

seq = CoefficientList[
AsymptoticDSolveValue[{eq == 0, y[0] == 1}, y, {x, 0, 30}],
x]

(*
{1, 1, 2, 5, 14, 46, 166, 652, 2780, 12644, 61136, 312676, 1680592, \
9467680, 55704104, 341185496, 2170853456, 14314313872, 97620050080, \
687418278544, 4989946902176, 37286121988256, 286432845428192, \
2259405263572480, 18280749571449664, 151561941235370176, \
1286402259593355776, 11168256342434121152, 99099358725069658880, \
898070590439513534464, 8306264068494786829696}
*)


And FindSequenceFunction returns the original recursion:

an = FindSequenceFunction[%]

(*
DifferenceRoot[
Function[{\[FormalY], \[FormalN]}, {-\[FormalN] (1 + \[FormalN]) \
\[FormalY][\[FormalN]] + (-2 - 2 \[FormalN]) \[FormalY][
1 + \[FormalN]] - 2 \[FormalY][2 + \[FormalN]] +
2 \[FormalY][3 + \[FormalN]] == 0, \[FormalY][1] ==
1, \[FormalY][2] == 1, \[FormalY][3] == 2}]]
*)

Table[an[n+1], {n, 0, 30}]

(*  same as above: 1, 1, 2, 5, 14...  *)


I guess I'm not sure what exactly is being sought, though.

• Nice soluton (+1). I tried AsymptoticDSolveValue[{x^5 y''[x]/2 + (2 x + 1) x^3 y'[ x] + (x^3 + x^2 + x - 1) y[x] + 1 == 0, y[0] == 1, y'[0] == 1}, y, {x, 0, 1}] instead, which returned unevaluated. Any idea why inclusion of y'[0] == 1 caused AsymptoticDSolveValue to fail? Oct 14, 2020 at 21:47
• @bbgodfrey I guess it has to do with $x=0$ being an irregular singular point. I'm not well-acquainted irregular singular points. But try a higher order, such as AsymptoticDSolveValue[{eq == 0, y[0] == 1, y'[0] == 1}, y, {x, 0, 2}]. (I don't actually recall removing y'[0] == 1, and it may have been carelessness on my part in copying code. But it worked and I didn't catch it.) Oct 15, 2020 at 1:49
• Amazing! {x, 0, 1} and x -> 0 don't work, but {x, 0, 2} and higher do. A bug? Oct 15, 2020 at 2:23
• @bbgodfrey I don't know (whether a bug). Possibly a heuristic about how far to develop the series to solve for the order 1 leads to some problem; however, SeriesTermGoal does not help. Oct 15, 2020 at 2:48
• @MichaelMykhaylov Did you "rationalize" the floating-point coefficients (i.e., change 0.5 to 1/2 and 1.5 to 3/2)? (bbgodfrey did this, too.) I gave the value of eq that I used in the first output line in my answer. You should use this and avoid floating point. Floating point cannot represent integers exactly above around 10^16, and Mathematica displays only the first six digits of FP numbers. You can see the complete number if you convert it to InputForm or FullForm`. Oct 15, 2020 at 13:35