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*}
