How do you calculate the integral of xy' = y + √(x^2+y^2) as a homogeneous differential equation?

First, let v = y/x

Now the equation is
y' = v + √(1+v^2)
Now, what about y'? That is dy/dx, but now we have a bunch of v stuff. Luckily,
y = vx
y' = xv' + v
Now we have
xv'+v = v + √(1+v^2)
xv' = √(1+v^2)
dv/√(1+v^2) = dx/x
arcsinh(v) = lnx+c
v = sinh(lnx+c)
Now, v=y/x, so
y = x sinh(lnx+c)
That's pretty simple, but we can get rid of the sinh at the cost of a messier answer.
Now, recall that sinh(z) = (e^z - e^-z)/2 so
y = x (e^(lnx+c) - e^-(lnx+c))/2
Now, since e^lnx = x, e^(lnx+c) = xe^c
y = x(xe^c - 1/(xe^c))/2
= x(x^2e^2c-1)/(2xe^c)
= ((xe^c)^2-1)/(2e^c)
I think I prefer the sinh form.