Let $$M(x, y) = 3x^2y + 2xy + y^3 \qquad \text{ and } \qquad N(x, y) = x^2 + y^2$$
Then, $$\frac{\partial}{\partial y}M(x, y) = 3x^2 + 2x + 3y^2 \qquad \text{ and } \qquad \frac{\partial}{\partial x}N(x, y) = 2x$$
We can see that this equation is not exact, however, note that $$\frac{M_y - N_x}{N} = \frac{3(x^2 + y^2)}{x^2 + y^2} = 3$$ is a function of $x$ only. Thus, there is an integrating factor $\mu(x)$ that satisfies the differential equation $$\frac{\mathrm{d}\mu}{\mathrm{d}x} = \left(\frac{M_y - N_x}{N}\right)\mu = 3\mu \qquad \implies \qquad \mu = e^{3x}$$
Multiplying our original equation by $\mu(x)$, we have $$(3x^2e^{3x}y + 2xe^{3x}y + e^{3x}y^3) + e^{3x}(x^2 + y^2)\frac{\mathrm{d}y}{\mathrm{d}x} = 0$$
We can see that this equation is exact, since $$\frac{\partial}{\partial y}(3x^2e^{3x}y + 2xe^{3x}y + e^{3x}y^3) = 3x^2e^{3x} + 2xe^{3x} + 3e^{3x}y^2 = \frac{\partial}{\partial x} e^{3x}(x^2 + y^2) $$
Thus, there exists a function $\psi(x, y)$ such that \begin{align*}\psi_x(x, y) &= 3x^2e^{3x}y + 2xe^{3x}y + e^{3x}y^3 \tag{1} \\\psi_y(x, y) &= e^{3x}(x^2 + y^2) \tag{2}\end{align*}
Integrating (1) with respect to $x$, we get \begin{align*}\psi(x, y) &= \int e^{3x}(3x^2y + 2xy + y^3) \\&= e^{3x}\left(x^2y + \frac{1}{3}y^3\right) + h(y) \tag*{by parts}\end{align*}
for some function $h$ of $y$. Next, differentiating with respect to $y$, and equating with (2), we get $$\psi_y(x, y) = e^{3x}(x^2 + y^2) + h'(y) = e^{3x}(x^2 + y^2)$$
Therefore, $$h'(y) = 0 \implies h \text{ is constant }$$ Taking $h(y) = 0$, we get $$\psi(x, y) = e^{3x}\left(x^2y + \frac{1}{3}y^3\right)$$
Thus, the solutions of the differential equation are given implicitly by $$e^{3x}\left(x^2y + \frac{1}{3}y^3\right) = C$$
