If you can assume that $\gamma$ is a small parameter, then write
$$x(t) = x_0(t) + \gamma x_1(t) + \gamma^2 x_2(t) + \cdots$$
$$y(t) = y_0(t) + \gamma y_1(t) + \gamma^2 y_2(t) + \cdots$$
i.e., assume such convergent series exist. A zeroth order solution is $y_0(t)=y(0)=y_0$, a constant, and the equation for $x_0(t)$ becomes
$$\frac{dx_0}{dt} = 2 (W+y_0) x_0-4 x_0^3$$
This equation is integrable:
$$\int \frac{dx_0}{2 (W+y_0) x_0-4 x_0^3} = t \implies \frac{1}{4 (W + y_0)} \log{\left (\frac{x_0(t)}{2 x_0(t)^2-W-y_0}\right)} = t+C$$
Solve for $x_0(t)$, then plug into the first order equation for $y_1(t)$:
$$\gamma \frac{dy_1}{dt} = \gamma(x_0(t)^2-y_0) \implies \frac{dy_1}{dt} = x_0(t)^2-y_0$$
Integrate with respect to $t$ to get $y_1(t)$, then plug into $x$ equation:
$$\frac{d}{dt} (x_0+\gamma x_1) = 2 (W + y_0+\gamma y_1) (x_0+\gamma x_1) - 4 (x_0+\gamma x_1)^3$$
Note that $(x_0+\gamma x_1)^3 = x_0^3 + 3 \gamma x_0^2 x_1 + O(\gamma^2)$. Coefficient of $\gamma^0$ is zero because of above equation. Equating coefficients of $\gamma^1$, we get
$$\frac{d x_1}{dt} = 2 (W+y_0) x_1 + x_0 y_1 - 12 x_0^2 x_1$$
This is an inhomogeneous 1st order equation for $x_1$, which may be solved using known techiques (e.g., integration factor).
At this point, you may repeat this process to get higher powers of $\gamma$. I do not have a proof that the resulting series converges.
Correction (after missing a sign:) As kobe pointed out, the original DE is $$
(x^2-y^2)y'-2xy=0,
$$
which as equation for a vector field reads
$$
(x^2-y^2)\,dy-2xy\,dx=0\iff Im(\bar z^2\,dz)=0\text{ with } z=x+iy.
$$
From the complex interpretation it is directly visible that this is not integrable, for that it would have to be an expression $Im(f(z)\,dz)$. But since $\bar z=|z|^2/z$, an integrating factor presents itself as $|z|^{-4}=(x^2+y^2)^{-2}$ which repairs this deficiency. Then
$$
\frac{(x^2-y^2)\,dy-2xy\,dx}{(x^2+y^2)^2}=Im(z^{-2}\,dz)=-d(Im(z^{-1}))
$$
which implies that all solution trajectories lie on the curves
$$
-Im(z^{-1})=\frac{y}{x^2+y^2}=C\quad(\in \Bbb R)
$$
which can be solved as
$$
y^2-2y(2C)^{-1}+(2C)^{-2}=(2C)^{-2}-x^2\\
y=(2C)^{-1}\pm\sqrt{(2C)^{-2}-x^2}
$$
giving the solutions $y\equiv0$ and $y=D\pm\sqrt{D^2-x^2}$ on the interval $x\in[-D,D]$.
Variant: If one reads too fast that the sign before $2xy$ were reversed, then
$$
(x^2-y^2)dy+2xydx=0\quad (\iff Im(z^2·dz)=0)
$$
is an exact differential equation with conserved quantity or first integral $x^2y-\frac13y^3$ ($=Im(\frac13 z^3)$) so that the solutions lie on the curves
$$
y(y^2-3x^2)=const.
$$
Best Answer
Hint. Assume that a polynomial of the form $y(x)=x^n+a_{n-1}x^{n-1}+\cdots+a_0$ satisfies the ODE, and pug it in the ODE. Then you shall find that $n=2$. Next, you shall find that $y=x^2-1/3$. Observe that if $y$ is a solution, then so is $cy$, for every $c\in \mathbb R$. Satisfaction of the condition $y(1)=2$, implies that the sought for solution is $$ y(x)=3x^2-1. $$