Exponential integral via a polynomial

$x^{2} y +xy^{2}+xy= x^{2}y +x\left(y^{2} +y\right) = \left( x\sqrt{y} + \frac{y^{2} +y }{2\sqrt{y}}\right)^{2} - \left(\frac{y^{2} +y }{2\sqrt{y}}\right)^{2}$

for simplicity let $\left(\frac{y^{2} +y }{2\sqrt{y}}\right)^{2} = m^{2}$ , a constant

we have $x^{2} y +xy^{2}+xy = \left( x\sqrt{y} + m\right)^{2} -m^{2}$

$\displaystyle \int_{-\infty}^{\infty}\left(\left( x^{2} y +xy^{2}+xy \right)e^{-\left(x^{2} y +xy^{2}+xy \right)} dx = \left(\left( x\sqrt{y} + m\right)^{2} -m^{2}\right)e^{-\left(\left( x\sqrt{y} + m\right)^{2} -m^{2}\right)} dx\right)$

$\displaystyle \int_{-\infty}^{\infty} \left(\left( x\sqrt{y} + m\right)^{2} -m^{2}\right)e^{-\left(\left( x\sqrt{y} + m\right)^{2} -m^{2}\right)} dx=\displaystyle \int_{-\infty}^{\infty}\left(x\sqrt{y} + m\right)^{2} e^{m^{2}} e^{-\left( x\sqrt{y} + m\right)^{2} }dx -\displaystyle \int_{-\infty}^{\infty}m^{2}e^{m^{2}} e^{-\left( x\sqrt{y} + m\right)^{2} }dx$

by transforming $x\sqrt{y} + m = u$ , we have

$\displaystyle \int_{-\infty}^{\infty}\frac{u^{2} e^{m^{2}} e^{-u^{2} }}{\sqrt{y}}du -\displaystyle \int_{-\infty}^{\infty}\frac{ m^{2}e^{m^{2}} e^{-u^{2} }}{\sqrt{y}}du$

we can solve each one separately:

$\displaystyle \int_{-\infty}^{\infty}\frac{ m^{2}e^{m^{2}} e^{-u^{2} }}{\sqrt{y}}du =\displaystyle \int_{-\infty}^{\infty} \frac{ m^{2}e^{m^{2}} e^{-u^{2} }}{\sqrt{y}} du$

This is a Guassian integral scaled by a factor of $\frac{ m^{2}e^{m^{2}} }{\sqrt{y}}$

$\displaystyle \int_{-\infty}^{\infty} \frac{e^{m^{2}} e^{-u^{2} }}{\sqrt{y}}du =\frac{m^{2}e^{m^{2}}}{\sqrt{y}} \times \sqrt{\pi}$

Similarly,

$\displaystyle \int_{-\infty}^{\infty}\frac{u^{2} e^{m^{2}} e^{-u^{2} }}{\sqrt{y}}du =\displaystyle \frac{ e^{m^{2}} }{\sqrt{y}} \int_{-\infty}^{\infty}u^{2}e^{-u^{2}} du$

this is $\frac{1}{2}$ the Gaussian integral ( you can verify it by using integration by parts), scaled by a factor of $\frac{ e^{m^{2}} }{\sqrt{y}}$

Therefore $\displaystyle \int_{-\infty}^{\infty}\frac{u^{2} e^{m^{2}} e^{-u^{2} }}{\sqrt{y}}du =\frac{\sqrt{\pi} e^{m^{2}} }{2\sqrt{y}}$

Now we have $\displaystyle \int_{-\infty}^{\infty}\frac{u^{2} e^{m^{2}} e^{-u^{2} }}{\sqrt{y}}du -\displaystyle \int_{-\infty}^{\infty}\frac{ m^{2}e^{m^{2}} e^{-u^{2} }}{\sqrt{y}}du =\frac{ \sqrt{\pi}e^{m^{2}} }{2\sqrt{y}}-\frac{\sqrt{\pi} m^{2}e^{m^{2}}}{\sqrt{y}} =0$

By rearranging we will get $m^{2}= \frac{1}{2}$

From our initial simplification $\left(\frac{y^{2} +y }{2\sqrt{y}}\right)^{2} = m^{2}$

$\left(\frac{y^{2} +y }{2\sqrt{y}}\right)^{2} = 1/2$

rearrange again to get $y^{2} +y = \sqrt{2y}$

Square the expression:

$y^{4} +2y^{3} +y^{2}=2y$

factor the $y$

$\displaystyle \mathbf{ y^{3} +2y^{2} +y=2}$

Did almost the same as Amal Hari's solution below. First complete the square in both the exponent and the factor: $x^2y+xy^2+xy = y\left(x^2 + x(y+1) \pm \frac{(y+1)^2}{4} \right) = y\left(x+\frac{y+1}{2}\right)^2 - \frac{y(y+1)^2}{4}=:yu^2-c$ Notice $4c$ is the value we need to find. Insert the simplifications and use $u$ -substitution with the shift above: $0 = \int_{\mathbb{R}}e^{-yu^2 + c}(yu^2-c)\:du = e^{c}\int_\mathbb{R} e^{-yu^2}(yu^2-c)\:du\qquad |:e^c\neq 0\quad (*)$ Use integration by parts to simplify the first integrand and get two "Gaussian Integrals": $\int_{\mathbb{R}}u^2e^{-yu^2}\:du=\left. -u\cdot\frac{e^{-yu^2}}{2y} \right|_{-\infty}^\infty +\frac{1}{2y}\int_\mathbb{R} 1\cdot e^{-yu^2}\:du=0+\frac{1}{2y}\cdot\sqrt{\frac{\pi}{y}}$ Insert the solutions of both "Gaussian Integrals" into our equation $(*)$ : $0=\sqrt{\frac{\pi}{y}}\left(\frac{\cancel{y}}{2\cancel{y}}-c\right)\quad\Rightarrow\quad y^3+2y^2+y=4c=\boxed{2}$