No Mathematica, você precisa colocá-lo da seguinte forma:
eqn = x*y'[x] - y[x] == Sqrt[1 + x^2]*y'[x] + Sqrt[1 + y[x]^2]
DSolve[eqn, y[x], x]
Com isso, obtenho:
{{y[x] -> InverseFunction[1/2 (ArcSinh[#1] + #1 (-#1 + Sqrt[1 + #1^2])) &][ 1/2 (-x (x + Sqrt[1 + x^2]) - ArcSinh[x]) + C[1]]}}
com o seguinte aviso:
Solve::tdep: The equations appear to involve the variables to be solved for in an essentially non-algebraic way.
Da ajuda, parece que o solucionador encontrou uma função transcendental. Isso é tanto quanto eu posso levar isso. Espero que isso ajude você a seguir em frente.