Solving the problem : $z z_x + z_y = 0, \quad z(x,0) = x^2$ Exercise :

For the problem :
  $$\begin{cases} zz_x + z_y = 0 \\ z(x,0) = x^2\end{cases}$$
  derive the solution :
  $$z(x,y) = \begin{cases} x^2, \quad y = 0\\ \frac{1+2xy - \sqrt{1+4xy}}{2y^2}, \quad y \neq 0 \; \text{and} \; 1+4xy >0 \end{cases}$$
  When do shocks develop ? Use the Taylor series for $\sqrt{1+\epsilon}$ about $\epsilon= 0$ to verify that $\lim_{y\to 0} z(x,y) = x^2$.

Attempt :
$$\frac{\mathrm{d}x}{z} = \frac{\mathrm{d}y}{1} = \frac{\mathrm{d}z}{0}$$
We yield the integral curves :
$$\frac{\mathrm{d}y}{1} = \frac{\mathrm{d}z}{0} \implies z_1 = z $$
$$\frac{\mathrm{d}x}{z} = \frac{\mathrm{d}y}{1} \implies z_2 = x -zy$$
Thus, the general solution will involve a $F \in C^1$ function, such that :
$$z(x,y) = F(x-zy)$$
For $y=0$ :
$$z(x,0) = F(x) \Rightarrow F(x) = x^2$$
How would one proceed now to find the second branch of the solution ?
(The shock-taylor part is easy)

The general solution is :
$$z(x,y) = F(x-zy)\qquad\text{OK}.$$
$F$ is an arbitrary function, to be determines according to the boundary condition.
Condition : $\quad z(x,0)=x^2=F(x-0y)$
So, the function $F$ is determined :
$$F(X)=X^2\qquad\text{any }X$$
We put this function into the general solution where $X=x-zy$.
$$z(x,y)=F(x-zy)=(x-zy)^2$$
$z=x^2-2xyz+y^2z^2$
$y^2z^2-(2xy+1)z+x^2=0$
Solving for $z$ leads to :
$$z=\frac{2xy+1\pm\sqrt{(2xy+1)^2-4x^2y^2}}{2y^2}$$
$$z=\frac{2xy+1\pm\sqrt{1+4xy}}{2y^2}$$
For $y\to 0$ :
Let $\quad 4xy=\epsilon>0 \quad$ because $1+4xy>0$.
$$z=\frac{\frac{\epsilon}{2}+1\pm\sqrt{1+\epsilon}}{2y^2}=\frac{\frac{\epsilon}{2}+1\pm\sqrt{1+\epsilon}}{2\left(\frac{\epsilon}{4x}\right)^2}=8x^2\left(\frac{\frac{\epsilon}{2}+1\pm\sqrt{1+\epsilon}}{\epsilon^2}\right)$$
$\sqrt{1+\epsilon}\simeq 1+\frac12\epsilon-\frac18\epsilon^2+...$
$$z\simeq 8x^2\left(\frac{\frac{\epsilon}{2}+1\pm\left(1+\frac12\epsilon-\frac18\epsilon^2+...\right)}{\epsilon^2}\right)$$
Case of sign $+$ :
$z\simeq 8x^2\left(\frac{2+\epsilon-\frac18\epsilon^2+...}{\epsilon^2}\right)\to\infty\quad$ when $\epsilon\to 0$. This case is rejected because $z\to x^2$.
Case of sign $-$ :
$z\simeq 8x^2\left(\frac{\frac18\epsilon^2+...}{\epsilon^2}\right)\to x^2\quad$when $\epsilon\to 0$. This agrees. Thus the final result is :
$$z=\frac{2xy+1-\sqrt{1+4xy}}{2y^2}$$
