Utilizar el cambio de variable $$ X=x^2,\quad Y=y^2,\quad U=u^2\\ P=\frac{\partial U}{\partial X},\quad Q=\frac{\partial U}{\partial Y} $$ Tenga en cuenta que $$ P=\frac{\partial U}{\partial u}\frac{\partial u}{\partial x}\frac{\partial x}{\partial X}=\frac{up}{x}\\ Q=\frac{\partial U}{\partial u}\frac{\partial u}{\partial y}\frac{\partial y}{\partial Y}=\frac{uq}{y} $$ Nuestro pde en términos de las nuevas variables es $$\begin{align} & P^2X+Q^2Y=X+Y\\\implies& X(P^2-1)+Y(Q^2-1)=0 \end{align}$$
Método general para so $$f(x,y,u,p,q)=f_1(x,p)+f_2(y,q)=0$$
En este caso las ecuaciones auxiliares de Charpit se convierten en $$ \frac{dp}{\frac{\partial f_1}{\partial x}}=\frac{dq}{\frac{\partial f_2}{\partial x}}=\frac{du}{-p\frac{\partial f_1}{\partial p}-q\frac{\partial f_2}{\partial q}}=\frac{dx}{-\frac{\partial f_1}{\partial p}}=\frac{dy}{-\frac{\partial f_2}{\partial q}} $$ Multiplicando por la 1ª y la 4ª relación obtenemos $$\begin{align} &\frac{\partial f_1}{\partial x}dx+\frac{\partial f_1}{\partial p}dp=0 \implies df_1=0\\ \therefore\;&f_1(x,p)=a \end{align}$$ donde $a$ es una constante. Y, $$f_2(y,q)=-f_1(x,p)=-a$$ Podemos escribir $p,q$ como $$p=F_1(x,a),\quad q=F_2(y,a)$$ Finalmente, integrando obtenemos la solución requerida $$u=\int F_1(x,a)\,dx+\int F_2(y,a)\,dy$$
Problema actual
$$\begin{align} F_1(X,a)=&\sqrt{1+\frac{a}{X}},\quad F_2(Y,a)=\sqrt{1-\frac{a}{Y}}\\ \implies U=&\int F_1(X,a)\,dX+\int F_2(Y,a)\,dY\\ \implies u^2=& 2\int\sqrt{x^2+a}\,dx+2\int\sqrt{y^2-a}\,dy\\ =& b+x\sqrt{x^2+a}+\sqrt a\ln|x+\sqrt{x^2+a}|+y\sqrt{y^2-a}-\sqrt a\ln|y+\sqrt{y^2-a}| \end{align}$$ Aquí hemos supuesto $a>0$ y $b$ es una constante arbitraria.