Assuming we're looking for a power series solution centered around [tex]x=0[/tex], take
[tex]y=\displaystyle\sum_{n\ge0}a_nx^n[/tex]
[tex]y'=\displaystyle\sum_{n\ge1}na_nx^{n-1}[/tex]
[tex]y''=\displaystyle\sum_{n\ge2}n(n-1)a_nx^{n-2}[/tex]
Substituting into the ODE yields
[tex]\displaystyle\sum_{n\ge2}n(n-1)a_nx^{n-2}+\sum_{n\ge1}na_nx^{n+1}-\sum_{n\ge0}a_nx^{n+1}=0[/tex]
The first series starts with a constant term; the second series starts at [tex]x^2[/tex]; the last starts at [tex]x^1[/tex]. So, extract the first two terms from the first series, and the first term from the last series so that each new series starts with a [tex]x^2[/tex] term. We have
[tex]\displaystyle\sum_{n\ge2}n(n-1)a_nx^{n-2}=2a_2+6a_3x+\sum_{n\ge4}n(n-1)a_nx^{n-2}[/tex]
[tex]\displaystyle\sum_{n\ge0}a_nx^{n+1}=a_0x+\sum_{n\ge1}a_nx^{n+1}[/tex]
Re-index the first sum to have it start at [tex]n=1[/tex] (to match the the other two sums):
[tex]\displaystyle\sum_{n\ge4}n(n-1)a_nx^{n-2}=\sum_{n\ge1}(n+3)(n+2)a_{n+3}x^{n+1}[/tex]
So now the ODE is
[tex]\displaystyle\left(2a_2+6a_3x+\sum_{n\ge1}(n+3)(n+2)a_{n+3}x^{n+1}\right)+\sum_{n\ge1}na_nx^{n+1}-\left(a_0x+\sum_{n\ge1}a_nx^{n+1}\right)=0[/tex]
Consolidate into one series starting [tex]n=1[/tex]:
[tex]\displaystyle2a_2+(6a_3-a_0)x+\sum_{n\ge1}\bigg[(n+3)(n+2)a_{n+3}+(n-1)a_n\bigg]x^{n+1}=0[/tex]
Suppose we're given initial conditions [tex]y(0)=a_0[/tex] and [tex]y'(0)=a_1[/tex] (which follow from setting [tex]x=0[/tex] in the power series representations for [tex]y[/tex] and [tex]y'[/tex], respectively). From the above equation it follows that
[tex]\begin{cases}2a_2=0\\6a_3-a_0=0\\(n+3)(n+2)a_{n+3}+(n-1)a_n=0&\text{for }n\ge2\end{cases}[/tex]
Let's first consider what happens when [tex]n=3k-2[/tex], i.e. [tex]n\in\{1,4,7,10,\ldots\}[/tex]. The recurrence relation tells us that
[tex]a_4=-\dfrac{1-1}{(1+3)(1+2)}a_1=0\implies a_7=0\implies a_{10}=0[/tex]
and so on, so that [tex]a_{3k-2}=0[/tex] except for when [tex]k=1[/tex].
Now let's consider [tex]n=3k-1[/tex], or [tex]n\in\{2,5,8,11,\ldots\}[/tex]. We know that [tex]a_2=0[/tex], and from the recurrence it follows that [tex]a_{3k-1}=0[/tex] for all [tex]k[/tex].
Finally, take [tex]n=3k[/tex], or [tex]n\in\{0,3,6,9,\ldots\}[/tex]. We have a solution for [tex]a_3[/tex] in terms of [tex]a_0[/tex], so the next few terms ([tex]k=2,3,4[/tex]) according to the recurrence would be
[tex]a_6=-\dfrac2{6\cdot5}a_3=-\dfrac2{6\cdot5\cdot3\cdot2}a_0=-\dfrac{a_0}{6\cdot3\cdot5}[/tex]
[tex]a_9=-\dfrac5{9\cdot8}a_6=\dfrac{a_0}{9\cdot6\cdot3\cdot8}[/tex]
[tex]a_{12}=-\dfrac8{12\cdot11}a_9=-\dfrac{a_0}{12\cdot9\cdot6\cdot3\cdot11}[/tex]
and so on. The reordering of the product in the denominator is intentionally done to make the pattern clearer. We can surmise the general pattern for [tex]n=3k[/tex] as
[tex]a_{3k}=\dfrac{(-1)^{k+1}a_0}{(3k\cdot(3k-3)\cdot(3k-2)\cdot\cdots\cdot6\cdot3\cdot(3k-1)}[/tex]
[tex]a_{3k}=\dfrac{(-1)^{k+1}a_0}{3^k(k\cdot(k-1)\cdot\cdots\cdot2\cdot1)\cdot(3k-1)}[/tex]
[tex]a_{3k}=\dfrac{(-1)^{k+1}a_0}{3^kk!(3k-1)}[/tex]
So the series solution to the ODE is given by
[tex]y=\displaystyle\sum_{n\ge0}a_nx^n[/tex]
[tex]y=a_1x+\displaystyle\sum_{k\ge0}\frac{(-1)^{k+1}a_0}{3^kk!(3k-1)}[/tex]
Attached is a plot of a numerical solution (blue) to the ODE with initial conditions sampled at [tex]a_0=y(0)=1[/tex] and [tex]a_1=y'(0)=2[/tex] overlaid with the series solution (orange) with [tex]n=3[/tex] and [tex]n=6[/tex]. (Note the rapid convergence.)
