I am late, but at some point, someone is going to want a coded solution that isn't rotated by 90Â° CCW.

Let $M(x,y) = x$ and let $N(x,y) = y + x^2y + y^3$. Then,

\begin{align*}

M_y(x,y) = 0, & N_x(x,y) = 2xy.

\end{align*}

So the equation is not exact, and we need an integrating factor $\mu(x,y) = \mu(y)$ to make it exact.

\begin{equation*}

\frac{d\mu(y)}{dy} = \frac{N_x - M_y}{M} \mu(y) = \frac{2xy}{x} \mu(y) = y\mu(y) \Longrightarrow \frac{d\mu(y)}{\mu(y)} = 2ydy \Longrightarrow \mu(y) = e^{y^2}

\end{equation*}

Now, we require a function $\psi(x,y)$ that satisfies $\psi_x = \mu M$ and $\psi_y = \mu N$.

\begin{equation*}

\psi_x(x,y) = xe^{y^2} \Longrightarrow \psi(x,y) = \frac{1}{2}x^2e^{y^2} + h(y)

\end{equation*}

We differentiate the result to get

\begin{equation*}

\psi_y(x,y) = x^2xe^{y^2} + h'(y) = e^{y^2}(y + x^2y + y^3) \Longleftrightarrow h'(y) = e^{y^2}(y^3 + y) \Longrightarrow h(y) = \int{(e^{y^2}(y^3 + y))}dy.

\end{equation*}

Let $u = y^2$ so that $du = 2ydy$. Then,

\begin{equation*}

h(y) = \int{(y^3e^{y^2})}dy + \int{(ye^{y^2})}dy = \frac{1}{2}\int{(ue^u)}du + \frac{1}{2}\int{e^u}du = \frac{1}{2}ue^u - \frac{1}{2}\int{e^u}du + \frac{1}{2}e^u = \frac{1}{2}ue^u - \frac{1}{2}e^u + \frac{1}{2}e^u = \frac{1}{2}y^2e^{y^2}.

\end{equation*}

Therefore, the solution is implicitly given by

\begin{equation*}

C = \frac{1}{2}e^{y^2}(x^2 + y^2).

\end{equation*}