\[ \left (x^2+1\right ) y(x)+y'(x)-y(x)^2-2 x=0 \] ✓ Mathematica : cpu = 0.0913241 (sec), leaf count = 49
\[\left \{\left \{y(x)\to \frac {e^{\frac {x^3}{3}+x}}{-\int _1^xe^{\frac {K[1]^3}{3}+K[1]}dK[1]+c_1}+x^2+1\right \}\right \}\] ✓ Maple : cpu = 0.09 (sec), leaf count = 34
\[y \relax (x ) = x^{2}+1+\frac {{\mathrm e}^{\frac {1}{3} x^{3}+x}}{c_{1}-\left (\int {\mathrm e}^{\frac {1}{3} x^{3}+x}d x \right )}\]
Hand solution
\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\relax (x) +Q\relax (x) y+R\relax (x) y^{2}\) where \(P\relax (x) =2x,Q\relax (x) =-\left (x^{2}+1\right ) ,R\relax (x) =1\). We can convert this to Bernoulli first order ODE in \(u\relax (x) \), 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\relax (x) \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\relax (x) =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\relax (x) & =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\relax (x) & =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*}