Theorem: Linear Second Order ODE/(x^2 - 1) y'' - 2 x y' + 2 y = (x^2 - 1)^2

The second order ODE:
:$(1): \quad \paren {x^2 - 1} y' ' - 2 x y' + 2 y = \paren {x^2 - 1}^2$
has the general solution:
:$y = C_1 x + C_2 \paren {x^2 + 1} + \dfrac {x^4} 6 - \dfrac {x^2} 2$

Proof:
$(1)$ can be manipulated into the form:
:$y' ' - \dfrac {2 x} {x^2 - 1} y' + \dfrac 2 {x^2 - 1} y = x^2 - 1$
It can be seen that this is a nonhomogeneous linear second order ODE in the form:
:$y' ' + \map P x y' + \map Q x y = \map R x$
where:
:$\map P x = -\dfrac {2 x} {x^2 - 1}$
:$\map Q x = \dfrac 2 {x^2 - 1}$
:$\map R x = x^2 - 1$
First we establish the solution of the corresponding homogeneous linear second order ODE:
:$\paren {x^2 - 1} y' ' - 2 x y' + 2 y = 0$
From Second Order ODE: $\paren {x^2 - 1} y' ' - 2 x y' + 2 y = 0$, this has the general solution:
:$y_g = C_1 x + C_2 \paren {x^2 + 1}$
It remains to find a particular solution $y_p$ to $(1)$.
Expressing $y_g$ in the form:
:$y_g = C_1 \map {y_1} x + C_2 \map {y_2} x$
we have:
{{begin-eqn}}
{{eqn | l = \map {y_1} x
      | r = x
      | c = 
}}
{{eqn | l = \map {y_2} x
      | r = x^2 + 1
      | c = 
}}
{{eqn | ll= \leadsto
      | l = \map { {y_1}'} x
      | r = 1
      | c = Power Rule for Derivatives
}}
{{eqn | l = \map { {y_2}'} x
      | r = 2 x
      | c = Power Rule for Derivatives
}}
{{end-eqn}}
By the Method of Variation of Parameters, we have that:
:$y_p = v_1 y_1 + v_2 y_2$
where:
{{begin-eqn}}
{{eqn | l = v_1
      | r = \int -\frac {y_2 \map R x} {\map W {y_1, y_2} } \rd x
      | c = 
}}
{{eqn | l = v_2
      | r = \int \frac {y_1 \map R x} {\map W {y_1, y_2} } \rd x
      | c = 
}}
{{end-eqn}}
where $\map W {y_1, y_2}$ is the Wronskian of $y_1$ and $y_2$.
We have that:
{{begin-eqn}}
{{eqn | l = \map W {y_1, y_2}
      | r = y_1 {y_2}' - y_2 {y_1}'
      | c = {{Defof|Wronskian}}
}}
{{eqn | r = x \paren {2 x} - \paren {x^2 + 1} \paren 1
      | c = 
}}
{{eqn | r = 2 x^2 - x^2 - 1
      | c = 
}}
{{eqn | r = x^2 - 1
      | c = 
}}
{{end-eqn}}
Hence:
{{begin-eqn}}
{{eqn | l = v_1
      | r = \int -\frac {y_2 \map R x} {\map W {y_1, y_2} } \rd x
      | c = 
}}
{{eqn | r = \int -\frac {\paren {x^2 + 1} \paren {x^2 - 1} } {x^2 - 1} \rd x
      | c = 
}}
{{eqn | r = -\int \paren {x^2 + 1} \rd x
      | c = 
}}
{{eqn | r = -\frac {x^3} 3 - x
      | c = Primitive of Power
}}
{{end-eqn}}
{{begin-eqn}}
{{eqn | l = v_2
      | r = \int \frac {y_1 \map R x} {\map W {y_1, y_2} } \rd x
      | c = 
}}
{{eqn | r = \int \frac {x \paren {x^2 - 1} } {x^2 - 1} \rd x
      | c = 
}}
{{eqn | r = \int x \rd x
      | c = 
}}
{{eqn | r = \frac {x^2} 2
      | c = Primitive of Power
}}
{{end-eqn}}
It follows that:
{{begin-eqn}}
{{eqn | l = y_p
      | r = \paren {-\frac {x^3} 3 - x} x + \paren {\frac {x^2} 2} \paren {x^2 + 1}
      | c = 
}}
{{eqn | r = -\frac {x^4} 3 - x^2 + \frac {x^4} 2 + \frac {x^2} 2
      | c = rearranging
}}
{{eqn | r = \frac {x^4} 6 - \frac {x^2} 2
      | c = rearranging
}}
{{end-eqn}}
So from General Solution of Linear 2nd Order ODE from Homogeneous 2nd Order ODE and Particular Solution:
:$y = y_g + y_p = C_1 x + C_2 \paren {x^2 + 1} + \dfrac {x^4} 6 - \dfrac {x^2} 2$
is the general solution to $(1)$.
{{qed}}