Another approach
With equations like these, it is often good to check for exactness
An equation in this form:
(x^2 + y^2) dx - 2xy dy = 0
....can be related to level surface \Phi(x,y) = const such that d\Phi = 0 = \Phi_x dx + \Phi_y dy
So in this case
\Phi_x =x^2 + y^2, \Phi_(xy) = 2y
\Phi_y =- 2xy, \Phi_(yx) = -2y
But because \Phi_(xy) ne \Phi_(yx) , the equation is not exact
We can then re-arrange the equation to say:
dy/dx = (x^2 +y^2)/(2xy) = x/(2y) + y/(2x)
This is a homogeneous equation, ie
dy/dx = f(x,y) = f(kx, ky) for any constant k
in this case we make a substitution v(x) = y/x, y = v x, so that y' = v' x + v and the new equation becomes separable
So
dy/dx = x/(2y) + y/(2x)
implies v'x + v = 1/(2v ) + v/2
v'x = (1-v^2)/(2v )
(2v)/(1-v^2) v' = 1/x
We integrate wrt x:
int( (2v)/(1-v^2) v' = 1/x ) dx
implies int (2v)/(1-v^2) dv =int 1/x dx
Or
-ln abs (1-v^2) = ln abs x + C \ \ \ [ = ln C abs x ] where, throughout, C is a generic constant
implies -1/ abs (1-v^2) = C abs x
For suitable v and x
1-v^2 = 1/ (C x)
1-(y/x)^2 = 1/ (C x)
y^2 =x^2 + C x
y = pm sqrt (x^2 + C x )