Solution (source code)

= Solution

Writing the two coefficients as $M,N$, their mixed <derivatives> satisfy $M_y=2xy=N_x$, so the smooth form is exact on $\mathbb R^2$. Integrating the first coefficient with respect to $x$ gives $f=2x^2+x^2y^2/2+h(y)$. Matching $f_y=y+x^2y$ gives $h'(y)=y$. One potential is therefore
$$
\boxed{f(x,y)=2x^2+\frac12(1+x^2)y^2}.
$$
The <integral curves> of $df=0$ are its <level curves>. Absorbing a factor of two into the constant gives
$$
4x^2+(1+x^2)y^2=C,\qquad
\boxed{y=\pm\sqrt{\frac{C-4x^2}{1+x^2}}}
$$
on real graph branches. The initial condition fixes $C=12$ and selects the positive branch:
$$
\boxed{y(x)=2\sqrt{\frac{3-x^2}{1+x^2}},\qquad-\sqrt3<x<\sqrt3}.
$$
This is the maximal smooth graph interval containing $x=1$; at either endpoint the <level curve> has a vertical tangent, so $y'$ cannot remain finite.