\[ \left (x^2+1\right ) y(x)+y'(x)-y(x)^2-2 x=0 \] ✓ Mathematica : cpu = 0.716725 (sec), leaf count = 48
\[\left \{\left \{y(x)\to \frac {e^{\frac {x^3}{3}+x}}{c_1-\int _1^x e^{\frac {K[1]^3}{3}+K[1]} \, dK[1]}+x^2+1\right \}\right \}\]
✓ Maple : cpu = 0.122 (sec), leaf count = 34
\[ \left \{ y \left ( x \right ) ={x}^{2}+1+{1{{\rm e}^{{\frac {{x}^{3}}{3}}+x}} \left ( {\it \_C1}-\int \!{{\rm e}^{{\frac {{x}^{3}}{3}}+x}}\,{\rm d}x \right ) ^{-1}} \right \} \]
\begin {align} \left ( x^{2}+1\right ) y+y^{\prime }-y^{2}-2x & =0\nonumber \\ y^{\prime } & =-\left ( x^{2}+1\right ) y+y^{2}+2x\tag {1} \end {align}
This is Riccati first order non-linear ODE of the form of the general form \(y^{\prime }=P\left ( x\right ) +Q\left ( x\right ) y+R\left ( x\right ) y^{2}\) where \(P\left ( x\right ) =2x,Q\left ( x\right ) =-\left ( x^{2}+1\right ) ,R\left ( x\right ) =1\). We can convert this to Bernoulli first order ODE in \(u\left ( x\right ) \), which is little easier to solve by using \(u=y-x^{2}-1\). The difference between Bernoulli and Riccati is that the term \(P\left ( x\right ) =0\) in Bernoulli. If \(P\left ( x\right ) \neq 0\) and \(R\left ( x\right ) \neq 0\) then it is called Riccati.
Using \(u=y-x^{2}-1\) gives\begin {align*} u^{\prime } & =y^{\prime }-2x\\ u^{\prime } & =\left [ -\left ( x^{2}+1\right ) y+y^{2}+2x\right ] -2x\\ & =-\left ( x^{2}+1\right ) \left ( u+x^{2}+1\right ) +\left ( u+x^{2}+1\right ) ^{2}\\ & =\left ( u+x^{2}+1\right ) \left [ \left ( u+x^{2}+1\right ) -\left ( x^{2}+1\right ) \right ] \\ & =\left ( u+x^{2}+1\right ) u\\ & =u^{2}+u\left ( 1+x^{2}\right ) \end {align*}
We see now this is Bernoulli since \(P\left ( x\right ) =0\). To solve Bernoulli we always start by dividing by \(u^{2}\) giving\[ \frac {u^{\prime }}{u^{2}}=1+\frac {1}{u}\left ( 1+x^{2}\right ) \] Next we let \(v=\frac {1}{u}\), hence \(v^{\prime }=-\frac {u^{\prime }}{u^{2}}\)therefore the above becomes\begin {align*} -v^{\prime } & =1+v\left ( 1+x^{2}\right ) \\ v^{\prime }+v\left ( 1+x^{2}\right ) & =-1 \end {align*}
Integrating factor is \(e^{\int \left ( 1+x^{2}\right ) dx}=e^{\left ( x+\frac {x^{3}}{2}\right ) }\), therefore\[ d\left ( e^{\left ( x+\frac {x^{3}}{2}\right ) }v\right ) =-e^{\left ( x+\frac {x^{3}}{2}\right ) }\] Integrating\begin {align*} e^{\left ( x+\frac {x^{3}}{2}\right ) }v & =-\int e^{\left ( x+\frac {x^{3}}{2}\right ) }dx+C\\ v\left ( x\right ) & =e^{-\left ( x+\frac {x^{3}}{2}\right ) }\left ( C-\int e^{\left ( x+\frac {x^{3}}{2}\right ) }dx\right ) \end {align*}
Therefore \[ u=\frac {1}{v}=\frac {e^{\left ( x+\frac {x^{3}}{2}\right ) }}{\left ( C-\int e^{\left ( x+\frac {x^{3}}{2}\right ) }dx\right ) }\] And since \(u=y-x^{2}-1\) then\begin {align*} y\left ( x\right ) & =u+1+x^{2}\\ & =\frac {e^{\left ( x+\frac {x^{3}}{2}\right ) }}{\left ( C-\int e^{\left ( x+\frac {x^{3}}{2}\right ) }dx\right ) }+1+x^{2} \end {align*}