Write the PDE as $F(x,y,u,u_{x},u_{y})=0$, where $F(x,y,z,p,q)=p^{2}+q^{2}-z^{2}$. The characteristic equations are
$$\begin{cases}\dot{x}(s)=F_{p}(x(s),y(s),z(s),p(s),q(s))=2p(s)\\ \dot{y}(s)=F_{q}(x(s),y(s),z(s),p(s),q(s))=2q(s)\\ \dot{z}(s)=p(s)F_{p}(x(s),y(s),z(s),p(s),q(s))+q(s)F_{q}(x(s),y(s),z(s),p(s),q(s))=2(p(s)^{2}+q(s)^{2})=2z(s)^{2}\\ \dot{p}(s)=-F_{x}(x(s),y(s),z(s),p(s),q(s))-p(s)F_{z}(x(s),y(s),z(s),p(s),q(s))=2p(s)z(s)\\ \dot{q}(s)=-F_{y}(x(s),y(s),z(s),p(s),q(s))-q(s)F_{z}(x(s),y(s),z(s),p(s),q(s))=2q(s)z(s)\end{cases}$$
Suppose that we have initial data $(x^{0},y^{0},z^{0},p^{0},q^{0})$. We first solve the ODE for $z(s)$, obtaining
$$z(s)=\dfrac{z^{0}}{1-2z^{0}s}$$
We substitute this result in to the ODEs for $p(s),q(s)$ to obtain
$$\begin{cases}\dot{p}(s)=2p(s)\dfrac{z^{0}}{1-2z^{0}s}, & p(s)=\dfrac{p^{0}}{1-2z^{0}s}\\ \dot{q}(s)=2q(s)\dfrac{z^{0}}{1-2z^{0}s}, & q(s)=\dfrac{q^{0}}{1-2z^{0}s} \end{cases}$$
Lastly, we substitute these results into the ODEs for $x(s),y(s)$ to obtain
$$\begin{cases}\dot{x}(s)=2\dfrac{p^{0}}{1-2z^{0}s}, & x(s)=-\dfrac{p^{0}}{z^{0}}\ln\left|1-2z^{0}s\right|+x^{0}\\ \dot{y}(s)=2\dfrac{q^{0}}{1-2z^{0}s}, & y(s)=-\dfrac{q^{0}}{z^{0}}\ln\left|1-2z^{0}s\right|+y^{0}\end{cases}$$
for $s$ in the appropriate domain.
Suppose our boundary condition is $u=1$ on $\Gamma=S^{1}$. Fix a point $(x,y)\in\mathbb{R}^{2}$, and let $s$ and $(x^{0},y^{0})=(\cos\theta^{0},\sin\theta^{0})$, for some $\theta^{0}\in[0,2\pi]$ be such that $(x,y)=(x(s),y(s))$. Clearly, $z^{0}=1$, so it remains to solve for $x^{0},y^{0},p^{0},q^{0},s$. Since $\frac{d}{d\theta}u(\cos\theta,\sin\theta)_{\theta^{0}}=Du(x^{0},y^{0})\cdot (-y^{0},x^{0})=0$, we see that $(p^{0},q^{0})=\lambda(x^{0},y^{0})$, where $\lambda$ is some nonzero scalar. Since ${p^{0}}^{2}+{q^{0}}^{2}={z^{0}}^{2}=1$, we see that $\left|\lambda\right|=1$.
$$x^{2}+y^{2}=\left({x^{0}}^{2}+{y^{0}}^{2}\right)\left[\lambda^{-1}(1-\lambda\ln\left|1-2s\right|)\right]^{2}\Rightarrow s=\dfrac{1-\exp[\lambda^{-1}(1-\sqrt{x^{2}+y^{2}})]}{2}$$
We conclude that
$$u(x,y)=u(x(s),y(s))=z(s)=\exp\left[-\lambda^{-1}(1-\sqrt{x^{2}+y^{2}})\right]=\exp\left[\pm(1-\sqrt{x^{2}+y^{2}})\right]$$
Suppose now our boundary condition is $u=1$ on $\Gamma=\mathbb{R}\times\left\{0\right\}$. Proceeding as above, we fix a point $(x,y)\in\mathbb{R}^{2}$, and let $s$ and $(x^{0},0)$ be such that $(x,y)=(x(s),y(s))$. As before, $z^{0}=1$. Since $u=1$ on $\Gamma$, we obtain that $p^{0}=u_{x}(x^{0},0)=0$. From the identity ${p^{0}}^{2}+{q^{0}}^{2}={z^{0}}^{2}$, we obtain that $\left|q^{0}\right|=1$. Substituting this information into the equations above, we see that
$$x=x^{0}, y=-q^{0}\ln\left|1-2s\right|\Rightarrow s=\dfrac{1-e^{\pm y}}{2}$$
Substituting this result into the equation for $z(s)$, we conclude that
$$u(x,y)=u(x(s),y(s))=z(s)=e^{\pm y}$$
Hint :
You got $p=\frac{c}{y^2}$. Put this in the given equation and find the value of $q$. Then use $\,dz=p\,dx+q\,dy $ and hence solve it.
Finally get , $$\,dz=\frac{c}{y^2}\,dx+\frac{c^2-2y^3z-2cxy}{2y^4}\,dy.$$
$$\implies2(y\,dz+z\,dy)=\frac{c^2}{y^3}\,dy+2c\frac{y\,dx-x\,dy}{y^2}$$
$$\implies 2\,d(zy)=\frac{c^2}{y^3}\,dy+2c\,d(x/y)$$
$$\implies2zy=-\frac{c^2}{2y^2}+2c.\frac{x}{y}+c_1$$
Best Answer
Obviously, the main difficulty comes from the boundary condition $u(x,0)=x^2+1$ because it involves to solve a polynomial equation leading to huge formula. Even if the solving is theoretically possible, one have to merely accept a result on implicit form.
To make more clear where the difficulty arises, we will avoid the profusion of symbols introduced into the usual method of characteristics, but in following the same approach in fact. $$u_x^2+u_y^2=u \quad;\quad u(x,0)=x^2+1$$ First change of function :$\quad u=v^2\quad\begin{cases}u_x=2vv_x \\ u_y=2vv_y\end{cases}\quad\to\quad v_x^2+v_y^2=\frac{1}{4}\quad;\quad v(x,0)=\sqrt{x^2+1}$
$$v_y=\sqrt{\frac{1}{4}-v_x^2} \quad\to\quad v_{xy}=\frac{v_xv_{xx}}{\sqrt{\frac{1}{4}-v_x^2}}$$ Second change of function : $\quad w=v_x\quad\to\quad w_{y}=\frac{ww_{x}}{\sqrt{\frac{1}{4}-w^2}}$ $$\sqrt{\frac{1}{4}-w^2}\:w_y-2w\:w_x=0 \quad;\quad w(x,0)=v_x(x,0)=\frac{x}{\sqrt{x^2+1}}$$
This is a first order PDE. The set of characteristic ODEs is : $\quad \frac{dy}{\sqrt{\frac{1}{4}-w^2}}=\frac{dx}{-2w}=\frac{dw}{0}$
First family of characteristic curves, coming from $dw=0 \quad\to\quad w=c_1$
Second family of characteristic curves, from $\frac{dy}{\sqrt{\frac{1}{4}-c_1^2}}=\frac{dx}{-2c_1} \quad\to\quad 2c_1y+\sqrt{\frac{1}{4}-c_1^2}x=c_2$
The general solution can be presented on various forms :$\begin{cases} \Phi\left(w\:,\: 2wy+\sqrt{\frac{1}{4}-w^2}\:x\right)=0\\ w=f\left(2wy+\sqrt{\frac{1}{4}-w^2}\:x\right)\\ 2wy+\sqrt{\frac{1}{4}-w^2}\:x=F(w)\end{cases}$
where $\Phi$ , $f$ , $F$ are any differentiable functions. Any one of these functions has to be determined according to the boundary condition.
$$w(x,0)=\frac{x}{\sqrt{x^2+1}}\quad\to\quad 2\frac{x}{\sqrt{x^2+1}}0+\sqrt{\frac{1}{4}-\left(\frac{x}{\sqrt{x^2+1}}\right)^2}\:x=F\left(\frac{x}{\sqrt{x^2+1}}\right)$$
Let $t=\frac{x}{\sqrt{x^2+1}} \quad\to\quad x=\frac{t}{\sqrt{1-t^2}} \quad\to\quad \sqrt{\frac{1}{4}-t^2}\:\frac{t}{\sqrt{1-t^2}} =F\left(t\right)$
Now, the function $F$ is determined. We put it into the above general solution :
$$2wy+\sqrt{\frac{1}{4}-w^2}\:x=\sqrt{\frac{1}{4}-w^2}\:\frac{w}{\sqrt{1-w^2}}$$
$$\frac{4}{\sqrt{1-4w^2}}\:y+\frac{x}{w}=\frac{1}{\sqrt{1-w^2}}$$
Solving this equation for $w$ leads to $\quad w(x,y)$
In fact, this is a four degree polynomial equation. One can solve it analytically, but this involves huge formulas. That is the hitch.
So, we let $w$ on the implicit form of the above equation and, from it, we consider that $w(x,y)$ is known.
$w=\frac{u_x}{2\sqrt{u}}\quad\to\quad \int \frac{u_x}{2\sqrt{u}}=\sqrt{u}=\int w(x,y)dx$ $$u(x,y)=\left(\int w(x,y)dx \right)^2$$