\[ a y(x) (y(x)-x)+y'(x)-1=0 \] ✓ Mathematica : cpu = 8.87254 (sec), leaf count = 120
\[\left \{\left \{y(x)\to \frac {c_1 \left (\sqrt {\frac {\pi }{2}} \sqrt {a} x e^{\frac {a x^2}{2}} \text {erf}\left (\frac {\sqrt {a} x}{\sqrt {2}}\right )+1\right )+a x e^{\frac {a x^2}{2}}}{a \left (\frac {\sqrt {\frac {\pi }{2}} c_1 e^{\frac {a x^2}{2}} \text {erf}\left (\frac {\sqrt {a} x}{\sqrt {2}}\right )}{\sqrt {a}}+e^{\frac {a x^2}{2}}\right )}\right \}\right \}\] ✓ Maple : cpu = 0.349 (sec), leaf count = 72
\[ \left \{ y \left ( x \right ) ={1 \left ( 2\,\sqrt {a}{{\rm e}^{-1/2\,a{x}^{2}}}+x \left ( \sqrt {\pi }{\it Erf} \left ( {\frac {\sqrt {2}x}{2}\sqrt {a}} \right ) \sqrt {2}a+2\,{a}^{3/2}{\it \_C1} \right ) \right ) \left ( \sqrt {\pi }{\it Erf} \left ( {\frac {\sqrt {2}x}{2}\sqrt {a}} \right ) \sqrt {2}a+2\,{a}^{3/2}{\it \_C1} \right ) ^{-1}} \right \} \]
\begin {align} y^{\prime }+ay\left ( y-x\right ) -1 & =0\nonumber \\ y^{\prime } & =1-\left ( ay^{2}-ayx\right ) \nonumber \\ & =1+ayx-ay^{2}\tag {1} \end {align}
This is Riccati first order non-linear ODE \(y^{\prime }=P\left ( x\right ) +A\left ( x\right ) y+R\left ( x\right ) y^{2}\) with \(P\left ( x\right ) =1,Q\left ( x\right ) =-ax,R\left ( x\right ) =-a\). We can convert Riccati to Bernoulli which is easier to solve using the substitution \(u=y-x\)\begin {align*} u^{\prime } & =y^{\prime }-1\\ & =\left ( 1+ayx-ay^{2}\right ) -1\\ & =\left ( 1+a\left ( u+x\right ) x-a\left ( u+x\right ) ^{2}\right ) -1\\ & =1+aux+ax^{2}-a\left ( u^{2}+x^{2}+2ux\right ) -1\\ & =1+aux+ax^{2}-au^{2}-ax^{2}-2aux-1\\ & =-aux-au^{2}\\ u^{\prime } & =-aux-au^{2} \end {align*}
This is of the form \(u^{\prime }=P\left ( x\right ) +Q\left ( x\right ) u+R\left ( x\right ) u^{2}\) and since \(P\left ( x\right ) =0\) then it is Bernoulli differential equation. (when \(P\left ( x\right ) \neq 0\) and \(R\left ( x\right ) \neq 0\) it is Riccati). To solve Bernoulli we always start by dividing by \(u^{2}\)\[ \frac {u^{\prime }}{u^{2}}=-\frac {ax}{u}-a \] Then we let \(\zeta =\frac {1}{u}\), hence \(\zeta ^{\prime }=-\frac {u^{\prime }}{u^{2}}\), therefore the above becomes\begin {align*} -\zeta ^{\prime } & =-ax\zeta -a\\ \zeta ^{\prime }-ax\zeta & =a \end {align*}
Integrating factor is \(e^{-\int axdx}=e^{-a\frac {x^{2}}{2}}\), hence \(d\left ( e^{-a\frac {x^{2}}{2}}\zeta \right ) =ae^{-a\frac {x^{2}}{2}}\). Integrating both sides gives\[ e^{-a\frac {x^{2}}{2}}\zeta =a\int e^{-a\frac {x^{2}}{2}}dx+C \] But \[ \int e^{-a\frac {x^{2}}{2}}dx=\sqrt {\frac {\pi }{2a}}\operatorname {erf}\left ( \sqrt {\frac {a}{2}}x\right ) \] Therefore\begin {align*} e^{-a\frac {x^{2}}{2}}\zeta & =a\sqrt {\frac {\pi }{2a}}\operatorname {erf}\left ( \sqrt {\frac {a}{2}}x\right ) +C\\ \zeta & =e^{a\frac {x^{2}}{2}}\left ( a\sqrt {\frac {\pi }{2a}}\operatorname {erf}\left ( \sqrt {\frac {a}{2}}x\right ) +C\right ) \end {align*}
Hence\begin {align*} u & =\frac {1}{\zeta }\\ & =e^{-a\frac {x^{2}}{2}}\left ( a\sqrt {\frac {\pi }{2a}}\operatorname {erf}\left ( \sqrt {\frac {a}{2}}x\right ) +C\right ) ^{-1} \end {align*}
Since \(u=y-x\) then\begin {align*} y & =u+x\\ & =e^{-a\frac {x^{2}}{2}}\left ( a\sqrt {\frac {\pi }{2a}}\operatorname {erf}\left ( \sqrt {\frac {a}{2}}x\right ) +C\right ) ^{-1}+x\\ & =\frac {e^{-a\frac {x^{2}}{2}}}{\sqrt {\frac {a\pi }{2}}\operatorname {erf}\left ( \sqrt {\frac {a}{2}}x\right ) +C}+x \end {align*}
Verification
eq:=diff(y(x),x)+a*y(x)*(y(x)-x)-1 = 0; sol:=exp(-a*x^2/2)/(sqrt(a*Pi/2)*erf(sqrt(a/2)*x)+_C1)+x; odetest(y(x)=sol,eq); 0