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

Question

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

Answer 1

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} *)

Answer 2

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.