Here's a nifty trick - you can transform any linear $n$th-order (in)homogeneous ordinary differential equation into a linear $1$st-order (in)homogeneous $n$-dimensional system of ODEs, which can be solved with the power of matrix exponentials. Behold:
Write $v=y'$, and reinterpret the differential equation as a first-order system:
$$\begin{pmatrix} y \\ v \end{pmatrix}' = \begin{pmatrix}0 & 1\\ -2(c+e^{A+Bx}) & 2a(b-x)\end{pmatrix} \begin{pmatrix} y \\ v \end{pmatrix}.$$
If we denote the column vector $(y,v)^T$ as $\vec{y}$, we can write this system as $\vec{y}'=P(x)\vec{y}$. The solution is actually the same as in the $1$-dimensional case, but you need the matrix exponential,
$$\vec{y} = \exp\left(\int_{x_0}^x P(\tau)d\tau\right) \vec{y}_0 $$
$$= \exp \begin{pmatrix}0 & x-x_0\\ \frac{2}{B}e^A(e^{Bx_0}-e^{Bx})+2c(x_0-x) & 2ab(x-x_0)+a(x_0^2-x^2)\end{pmatrix} \vec{y}_0. $$
Before we actually attempt to synthesize the incoming monstrosity, let's prepare ourselves by writing this matrix using variables: $\begin{pmatrix} 0 & \alpha \\ \beta & \gamma\end{pmatrix}$. Because, let's face it, the full analytical solution is going to be insanely complicated. Also, we are going to do it using the Putzer Algorithm (just keep in mind that we need to distinguish between $x$ and $t$; while using this method we must fix $x$). The eigenvalues are $\lambda_{1,2} = (\gamma\pm\sqrt{\gamma^2+4\alpha\beta})/2$. Solving for $p_1$ and $p_2$ using generic methods, and plugging them into the formula with $t=1$, we get
$$\vec{y}=\left[ \begin{pmatrix} e^{\lambda_1} & 0 \\ 0 & e^{\lambda_2}\end{pmatrix} +\frac{e^{\lambda_1}-e^{\lambda_2}}{\lambda_1-\lambda_2} \begin{pmatrix} -\lambda_1 & \alpha \\ \beta & \gamma-\lambda_1 \end{pmatrix}\right]\vec{y}_0.$$
Just remember the $y$ you want is the first component of $\vec{y}$, and that $\lambda_{1,2}$ are functions of $\alpha,\beta,\gamma$, and that $\alpha,\beta,\gamma$ are all functions of $x$, so this technically contains your full analytical solution, just with a lot of substitutions. Have fun with your research OP. :)
$$y(1+xy)dx + x(1+xy+x^2y^2)dy = 0\tag1$$
Comparing $(1)$ with $~Mdx+Ndy=0~$, we have
$$M=y(1+xy)\qquad\text{and}\qquad N=x(1+xy+x^2y^2)$$
showing that $(1)$ is of the form $$f_1(xy)~y~dx+f_2(xy)~x~dy=0~.$$
Again $~Mx-Ny=xy(1+xy)-xy(1+xy+x^2y^2)=-x^3y^3\ne 0~$.
Hence the integrating factor (I.F.) is $~\dfrac{1}{Mx-Ny}=-\dfrac{1}{x^3y^3}~$.
Multiplying both side of $(1)$ by I.F. we have
$$\dfrac{1}{x^3y^2}(1+xy)dx + \dfrac{1}{x^2y^3}(1+xy+x^2y^2)dy = 0$$
$$\left[\dfrac{1}{x^3y^2}+\dfrac{1}{x^2y}\right]dx + \left[\dfrac{1}{x^2y^3}+\dfrac{1}{xy^2}+\dfrac{1}{y}\right]dy = 0$$
$$d\left(-\dfrac{1}{2x^{2}y^{2}}-\dfrac1{xy}+\ln y\right) = 0$$
Integrating we have
$$-\dfrac{1}{2x^{2}y^{2}}-\dfrac1{xy}+\ln y=c$$where $~c~$ is integrating constant.
Best Answer
The equation can be rewritten as $(y')^2+yy'-x^2-x=0.$ One can solve $y'$ in terms of $y$ and $x.$ Notice that the above equation is equivalent to $(2y')^2+2y(2y')-4x^2-4x=0=(2y')^2+2y(2y')+y^2-(y^2+4x^2+4x)=(2y'+y)^2-(y^2+4x^2+4x),$ hence $(2y'+y)^2=y^2+4x^2+4x,$ implying $y^2\geq-(4x^2+4x).$ Notice that $4x^2+4x=4x^2+4x+1-1=(2x+1)^2-1,$ so $y^2\geq1-(2x+1)^2.$ This allows us to say that $2y'+y=\pm\sqrt{y^2+(2x+1)^2-1}.$