\[ y''(x)-\left (x^2+1\right ) y(x)=0 \] ✓ Mathematica : cpu = 0.00783961 (sec), leaf count = 33
\[\left \{\left \{y(x)\to c_1 D_{-1}\left (\sqrt {2} x\right )+c_2 D_0\left (i \sqrt {2} x\right )\right \}\right \}\]
✓ Maple : cpu = 0.094 (sec), leaf count = 17
\[ \left \{ y \left ( x \right ) ={{\rm e}^{{\frac {{x}^{2}}{2}}}} \left ( {\it Erf} \left ( x \right ) {\it \_C2}+{\it \_C1} \right ) \right \} \]
\begin {equation} y^{\prime \prime }-\left ( x^{2}+1\right ) y=0 \tag {1} \end {equation}
Second order with varying coefficient. Using power series, let \(y=\sum _{n=0}^{\infty }c_{n}x^{n}\), hence\begin {align*} y^{\prime } & =\sum _{n=0}^{\infty }nc_{n}x^{n-1}=\sum _{n=1}^{\infty }nc_{n}x^{n-1}=\sum _{n=0}^{\infty }\left ( n+1\right ) c_{n+1}x^{n}\\ y^{\prime \prime } & =\sum _{n=0}^{\infty }n\left ( n+1\right ) c_{n+1}x^{n-1}=\sum _{n=1}^{\infty }n\left ( n+1\right ) c_{n+1}x^{n-1}=\sum _{n=0}^{\infty }\left ( n+1\right ) \left ( n+2\right ) c_{n+2}x^{n} \end {align*}
Substituting back in the original ODE gives\begin {align*} \sum _{n=0}^{\infty }\left ( n+1\right ) \left ( n+2\right ) c_{n+2}x^{n}-\left ( x^{2}+1\right ) \sum _{n=0}^{\infty }c_{n}x^{n} & =0\\ \sum _{n=0}^{\infty }\left ( n+1\right ) \left ( n+2\right ) c_{n+2}x^{n}-x^{2}\sum _{n=0}^{\infty }c_{n}x^{n}-\sum _{n=0}^{\infty }c_{n}x^{n} & =0\\ \sum _{n=0}^{\infty }\left ( n+1\right ) \left ( n+2\right ) c_{n+2}x^{n}-\sum _{n=0}^{\infty }c_{n}x^{n+2}-\sum _{n=0}^{\infty }c_{n}x^{n} & =0\\ \sum _{n=0}^{\infty }\left ( n+1\right ) \left ( n+2\right ) c_{n+2}x^{n}-\sum _{n=2}^{\infty }c_{n-2}x^{n}-\sum _{n=0}^{\infty }c_{n}x^{n} & =0 \end {align*}
For \(n=0\)\begin {align*} \left ( n+1\right ) \left ( n+2\right ) c_{n+2}-c_{n} & =0\\ 2c_{2}-c_{0} & =0\\ c_{2} & =\frac {c_{0}}{2} \end {align*}
For \(n=1\)\begin {align*} \left ( n+1\right ) \left ( n+2\right ) c_{n+2}-c_{n} & =0\\ \left ( 2\right ) \left ( 3\right ) c_{3}-c_{1} & =0\\ c_{3} & =\frac {c_{1}}{6} \end {align*}
For \(n\geq 2\)\begin {align*} \left ( n+1\right ) \left ( n+2\right ) c_{n+2}-c_{n-2}-c_{n} & =0\\ c_{n+2} & =\frac {c_{n-2}+c_{n}}{\left ( n+1\right ) \left ( n+2\right ) } \end {align*}
Hence for \(n=2\)\[ c_{4}=\frac {c_{0}+c_{2}}{\left ( 3\right ) \left ( 4\right ) }=\frac {c_{0}+\frac {c_{0}}{2}}{\left ( 3\right ) \left ( 4\right ) }=\frac {2c_{0}+c_{0}}{\left ( 2\right ) \left ( 3\right ) \left ( 4\right ) }=c_{0}\frac {3}{\left ( 2\right ) \left ( 3\right ) \left ( 4\right ) }\] For \(n=3\)\[ c_{5}=\frac {c_{1}+c_{3}}{\left ( 4\right ) \left ( 5\right ) }=\frac {c_{1}+\frac {c_{1}}{6}}{\left ( 4\right ) \left ( 5\right ) }=\frac {6c_{1}+c_{1}}{\left ( 4\right ) \left ( 5\right ) \left ( 6\right ) }=c_{1}\frac {7}{\left ( 4\right ) \left ( 5\right ) \left ( 6\right ) }\] For \(n=4\)\[ c_{6}=\frac {c_{2}+c_{4}}{\left ( 5\right ) \left ( 6\right ) }=\frac {\frac {c_{0}}{2}+c_{0}\frac {3}{\left ( 2\right ) \left ( 3\right ) \left ( 4\right ) }}{\left ( 5\right ) \left ( 6\right ) }=\frac {c_{0}\left ( 3\right ) \left ( 4\right ) +3c_{0}}{\left ( 2\right ) \left ( 3\right ) \left ( 4\right ) \left ( 5\right ) \left ( 6\right ) }=c_{0}\frac {15}{\left ( 2\right ) \left ( 3\right ) \left ( 4\right ) \left ( 5\right ) \left ( 6\right ) }\] For \(n=5\)\begin {align*} c_{7} & =\frac {c_{3}+c_{5}}{\left ( 6\right ) \left ( 7\right ) }\\ & =\frac {\frac {c_{1}}{6}+c_{1}\frac {7}{\left ( 4\right ) \left ( 5\right ) \left ( 6\right ) }}{\left ( 6\right ) \left ( 7\right ) }=\frac {3}{560}c_{1} \end {align*}
And so on. Hence the series is\begin {align*} y & =\sum _{n=0}^{\infty }c_{n}x^{n}\\ & =c_{0}+c_{1}x+c_{2}x^{2}+c_{3}x^{3}+\cdots \\ & =c_{0}+c_{1}x+\frac {c_{0}}{2}x^{2}+\frac {c_{1}}{6}x^{3}+c_{0}\frac {3}{\left ( 2\right ) \left ( 3\right ) \left ( 4\right ) }x^{4}+c_{1}\frac {7}{\left ( 4\right ) \left ( 5\right ) \left ( 6\right ) }x^{5}+c_{0}\frac {15}{\left ( 2\right ) \left ( 3\right ) \left ( 4\right ) \left ( 5\right ) \left ( 6\right ) }x^{6}+c_{1}\frac {3}{560}x^{7}+\cdots \\ & =c_{0}+c_{1}x+\frac {c_{0}}{2}x^{2}+\frac {c_{1}}{6}x^{3}+c_{0}\frac {1}{8}x^{4}+c_{1}\frac {7}{120}x^{5}+c_{0}\frac {1}{48}x^{6}+c_{1}\frac {3}{560}x^{7}+\cdots \\ & =c_{0}\left ( 1+\frac {x^{2}}{2}+\frac {1}{8}x^{4}+\frac {1}{48}x^{6}+\cdots \right ) +c_{1}\left ( x+\frac {1}{6}x^{3}+\frac {7}{120}x^{5}+\frac {3}{560}x^{7}+\cdots \right ) \end {align*}
Now the power series for \(e^{\frac {x^{2}}{2}}=1+\frac {x^{2}}{2}+\frac {x^{4}}{8}+\frac {x^{6}}{48}+\cdots \), so we can convert the first term above (the expression for \(c_{0}\) to be \(e^{\frac {x^{2}}{2}}\). Hence\[ c_{0}\left ( 1+\frac {x^{2}}{2}+\frac {1}{8}x^{4}+\frac {1}{48}x^{6}+\cdots \right ) =c_{0}e^{\frac {x^{2}}{2}}\] So now we have to work on the second term (the expression for \(c_{1}\))\[ c_{1}\left ( x+\frac {1}{6}x^{3}+\frac {7}{120}x^{5}+\frac {3}{560}x^{7}+\cdots \right ) =? \]
Recall that series for error function is\[ \operatorname {erf}\left ( x\right ) =\frac {2}{\sqrt {\pi }}\left ( x-\frac {x^{3}}{3}+\frac {x^{5}}{10}-\frac {x^{7}}{48}+\cdots \right ) \] Multiplying \(e^{\frac {x^{2}}{2}}\) by \(\operatorname {erf}\left ( x\right ) \) gives\begin {align*} e^{\frac {x^{2}}{2}}\operatorname {erf}\left ( x\right ) & =\frac {2}{\sqrt {\pi }}\left ( 1+\frac {x^{2}}{2}+\frac {x^{4}}{\left ( 2\right ) \left ( 4\right ) }+\frac {x^{6}}{\left ( 2\right ) \left ( 4\right ) \left ( 6\right ) }+\cdots \right ) \left ( x-\frac {x^{3}}{3}+\frac {x^{5}}{10}-\frac {x^{7}}{48}+\cdots \right ) \\ & =\frac {2}{\sqrt {\pi }}\left ( x-\frac {x^{3}}{3}+\frac {x^{5}}{10}-\frac {x^{7}}{48}+\cdots +\frac {x^{3}}{2}-\frac {x^{5}}{6}+\frac {x^{7}}{20}-\frac {x^{15}}{96}+\cdots \right ) \\ & =\frac {2}{\sqrt {\pi }}\left ( x+\frac {x^{3}}{6}+\frac {7x^{5}}{120}+\frac {3}{560}x^{7}+\cdots \right ) \end {align*}
Comparing the above to the term next to \(c_{1}\) above, we see they are the same with a multiplier \(\frac {2}{\sqrt {\pi }}\), which can be absorbed into the constant \(c_{1}\), Hence\begin {align*} y & =c_{0}\left ( 1+\frac {x^{2}}{2}+\frac {1}{8}x^{4}+\frac {1}{48}x^{6}+\cdots \right ) +c_{1}\left ( x+\frac {1}{6}x^{3}+\frac {7}{120}x^{5}+\frac {3}{560}x^{7}+\cdots \right ) \\ & =c_{0}e^{\frac {x^{2}}{2}}+c_{1}\left ( e^{\frac {x^{2}}{2}}\operatorname {erf}\left ( x\right ) \right ) \end {align*}
Hence final solution is\[ y=c_{0}e^{\frac {x^{2}}{2}}+c_{2}\left ( e^{\frac {x^{2}}{2}}\operatorname {erf}\left ( x\right ) \right ) \]
Verification
restart; ode:=diff(diff(y(x),x),x)-(x^2+1)*y(x)=0; y0:=_C1*exp(x^2/2)+_C2*exp(x^2/2)*erf(x); odetest(y(x)=y0,ode); 0