Let y=vx, where v is a function of x, which means dy/dx= v +x (dv)/dx. The given DE , thus transforms to v+x(dv)/dx =(v^3 +3v)/(1+3v^2)
x(dv)/dx = (v^3+3v)/(1+3v^2) -v =(2v-2v^3)/(1+3v^2)
(1+3v^2)/(v(1-v^2))dv= 2(1/x) dx. Partial fractions on the left side would be
(1/v +(4v)/(1-v^2))dv= 2 (1/x) dx. Now integrate on both sides to get
ln v -2ln(1-v^2) = 2 lnx +C. On simplification it becomes,
ln (v/(1-v^2)^2) = ln x^2 +C
v/(1-v^2)^2 = c_1x^2
(yx^3)/(x^2-y^2)^2= c_1x^2
(xy)/(x^2-y^2)^2 =c_1