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.