Let y=∑aᵣxʳ, then:

y'=∑raᵣxʳ⁻¹ and y''=∑r(r-1)aᵣxʳ⁻², so the DE can be written:

(x²+1)∑r(r-1)aᵣxʳ⁻²+x∑raᵣxʳ⁻¹-∑aᵣxʳ=0.

Since we have a 2nd order DE we can prepare two arbitrary constants A and B, which can be allocated as required.

So,

∑r(r-1)aᵣxʳ+∑raᵣxʳ-∑aᵣxʳ+∑r(r-1)aᵣxʳ⁻²=0.

We can collect xʳ terms together:

∑(r²-r+r-1)aᵣxʳ+∑r(r-1)aᵣxʳ⁻²=0,

∑(r²-1)aᵣxʳ+∑r(r-1)aᵣxʳ⁻²=0.

Then, by replacing r with r+2 in the second summation, we get:

∑(r²-1)aᵣxʳ+∑(r+2)(r+1)aᵣ₊₂xʳ=0.

Therefore:

∑[(r²-1)aᵣ+(r+2)(r+1)aᵣ₊₂]xʳ=0,

∑(r+1)[(r-1)aᵣ+(r+2)aᵣ₊₂]xʳ=0.

The coefficients of xʳ must all be zero (because this sum must be zero for all x), so we can find what aᵣ must be to satisfy this requirement.

Let r=0 (for the constant term):

-a₀+2a₂=0, a₂=a₀/2. We can use constant A instead of a₀, so a₂=A/2.

When r=1, 6a₃=0, and a₃=0. We can’t determine a₁ so we’ll use constant B in place of a₁.

We have an iterative equation for the a’s:

(r-1)aᵣ+(r+2)aᵣ₊₂=0.

This can be rewritten: (r-3)aᵣ₋₂+raᵣ=0,

aᵣ=(3-r)aᵣ₋₂/r.

If r=2, a₂=a₀/2=A/2, which we established earlier.

If r=4, a₄=-a₂/4=-A/8.

When r=3, a₃=0, so a₅=a₇=...=0.

When r=6, a₆=-3a₄/6=A/16,

r₈=-5a₆/8=-5A/128,

r₁₀=-7a₈/10=7A/256, and so on.

Therefore y=A+Bx+Ax²/2-Ax⁴/8+Ax⁶/16-5Ax⁸/128+7Ax¹⁰/256...,

y=Bx+A(1+x²/2-x⁴/8+x⁶/16-5x⁸/128+7x¹⁰/256...)

CHECK

y'=B+Ax-Ax³/2+3Ax⁵/8-5Ax⁷/16+35x⁹/128...

y''=A-3Ax²/2+15Ax⁴/8-35Ax⁶/16+315x⁸/128...

Now let’s assemble the DE:

(x²+1)y''=Ax²-3Ax⁴/2+15Ax⁶/8-35Ax⁸/16+315x¹⁰/128...

+A-3Ax²/2+15Ax⁴/8-35Ax⁶/16+315x⁸/128...

xy'=Bx+Ax²-Ax⁴/2+3Ax⁶/8-5Ax⁸/16+35x¹⁰/128...

-y=-A-Bx-Ax²/2+Ax⁴/8-Ax⁶/16+5Ax⁸/128-7Ax¹⁰/256...

Adding:

(A-A)+(B-B)x+(A-3A/2+A-A/2)x²+(-3A/2+15A/8-A/2+A/8)x⁴+...

You can see that each coefficient is zero.