If complete integral of differential equation $$ x (p^2 +q^2) = zp $$ ( p is partial derivative of z with respect to x and q is partial derivative of z with respect to y )
Passing through $x=0$ and $z^2 =4y$ then envelope of this family passing through $x=1$ ,$y=1 $ has
- 1) $z= -2 $
- 2) $z=2$
- 3) $z= √(2+2√2)$
- 4) $z= -√(2+2√2)$
As per https://en.wikipedia.org/wiki/Method_of_characteristics#Fully_nonlinear_case, a generalized form of characteristics is obtained by solving the ODE system (Lagrange-Charpit equations) \begin{align} \frac{dx}{2px-z} =\frac{dy}{2xq} =\frac{dz}{x(p^2+q^2)} =\frac{dz}{pz} =\frac{dp}{-q^2} =\frac{dq}{pq} \end{align} The last two ratios give $p^2+q^2=a^2=const.$, thus $a^2x=pz$.
The third last and last ratios give $\dfrac{dz}z=\dfrac{dq}q$ thus $$z=cq\tag{$1z$}$$ which gives $$x=\dfrac{c}{a^2}pq.\tag{$1x$}$$
Then also $\dfrac{dy}{2x}=\dfrac{dq}{p}=\dfrac{cq\,dq}{x}$ which gives $$y=\dfrac{c}{a^2}q^2+d.\tag{$1y$} $$
At this point the characteristic curve is represented as the image of the circle $p^2+q^2=a^2$. Elimination of the partial derivatives $p,q$ implies $$ (ax)^2+(a(y-d))^2=z^2\tag{2} $$ as the relation between $x,y,z$ along a characteristic curve.
The initial curve $x_0=0$, $y_0=\frac14z_0^2$ resp $s\mapsto(0,\frac14s^2,s)$ is part of the solution surface, thus $$ \dot z_0(s)=p_0(s)\dot x_0(s)+q_0(s)\dot y_0(s) \\ \implies 1=p_0·0+q_0·\frac12z_0~~\text{ or }~~q_0z_0=2\tag{3} $$ and from the original equation $$0·(p_0^2+q_0^2)=p_0z_0$$ which requires $p_0=0$ in the general case, thus $q_0^2=a^2$ and $ca^2=cq_0^2=q_0z_0=2$.
From $x=\dfrac{c}{a^2}pq=\dfrac12c^2pq$ we get $$ 4x^2=z^2(z_0^2-z^2)\tag{$4x$} $$ Also $y_0=c+d$ so that $$ y=\frac{z^2}{2}+y_0-\frac{z_0^2}2=\frac{z^2}{2}-\frac{z_0^2}4\tag{$4y$} $$
Thus to find the values on the characteristic curve through $x=1=y$ one needs to solve the system $$ 4=z_0^2z^2-z^4,\\ 4=2z^2-z_0^2. $$ for $z$.