Ok, I figured this out. It is clear that $y/z=\mathrm{const}=\psi$ will work. For $\phi$ we can substitute $z=y/c$. We end up with an exact equation:
$$2xy dx + (x^2+3y^2+\dfrac{3}{c^2}y^2) dy = \dfrac{d}{dy}(\phi(x(y),y)) = 0$$
with $$\phi=x^2y+y^3+\dfrac{yy^2}{c^2} = x^y+y^3+z^2y.$$
Clearly $\phi$ is a constant since its derivative is 0.
Thus: $F\left(\dfrac{u}{y},x^2y+y^3+u^2y\right)=0$ is the general solution for an arbitrary $F$ with continuous first derivatives. Using the implicit function theorem this is the same as:
$$u=yf(x^2y+y^3+u^2y).$$
Modifying the problem. Rather than consider the PDE $u_{xt}+u\ u_{xx}+\frac{1}{2}u_{x}^2=0$ with initial condition $u(x,0)=u_0(x)$ as asked above, I will consider the following variant.
$$
\text{Solve }u_{xy}+u\ u_{xx} + u_x^2=0\text{ subject to }u(x,0)=f(x).\qquad(\star)
$$
There are three differences between this question and that which was asked originally.
- The coefficient of $u_x^2$ has changed from $\frac{1}{2}$ to $1$.
- The variable $t$ has been renamed to $y$.
- The initial function $u_0(x)$ has been renamed to $f(x)$.
Only (1) represents a significant modification of the problem. It makes the solution more tractable and enables it to be found using an elementary application of the method of characteristics. For these reasons, it is conceivable that this was the intended question.
Note: I will not delve into regularity of the solutions in this answer.
Reduction to a first order quasilinear PDE. Write the equation as
$$
\frac{\partial}{\partial x}\left(u_y+u\ u_x\right)=0.
$$
Thus $(\star)$ is equivalent to
$$
u_y+u\ u_x=g(y),\qquad u(x,0)=f(x),\qquad (\star\star)
$$
where $g(y)$ is an arbitrary function of $y$ (with sufficient regularity).
Method of characteristics.
Perhaps the simplest formulation of the method of characteristics is for quasilinear first order PDEs. These are PDEs of the form
$$a(x,y,u)u_x+b(x,y,u)u_y=c(x,y,u).$$
To solve this equation, one regards the solution as a surface $z=u(x,y)$ in $xyz$-space. Let $s$ parametrize the initial curve $\bigl(s,0,f(s)\bigr)$ and let $t$ be a second parameter, which can be thought of as the distance flowed along a characteristic curve emanating from $\bigl(s,0,f(s)\bigr)$.
The characteristic equations are then
$$
\frac{dx}{dt}=a(x,y,z),\quad \frac{dy}{dt}=b(x,y,z),\quad \frac{dz}{dt}=c(x,y,z).
$$
Returning to our equation $(\star\star)$, this reduces to $a(x,y,u)=u$ and $b(x,y,u)=1$ and $c(x,y,u)=g(y)$.
Thus
$$
\frac{dx}{dt}=z,\quad \frac{dy}{dt}=1,\quad \frac{dz}{dt}=g(y)
$$
with initial conditions $x(0)=s$ and $y(0)=0$ and $z(0)=f(s)$.
The solution to this system is
$$
x=s+zt,\quad y=t,\quad z=f(s)+h(t),
$$
where $h(t)$ is the antiderivative of $g(t)$ satisfying $h(0)=0$. Since $g$ was arbitrary, so is $h$ given $h(0)=0$.
The solution. Now we eliminate all occurrences of $t$ by replacing them with $y$, then eliminate $s$ by writing $s=x-zy$. Finally, replace $z$ with $u$ to obtain the implicit equation
$$
\boxed{u=f(x-uy)+h(y)},
$$
where $h(y)$ is any sufficiently regular function satisfying $h(0)=0$. This is an implicit equation for the general solution of $(\star)$.
TL;DR. Change the $\frac{1}{2}$ in the original question to $1$ to obtain a PDE solvable by the method of characteristics.
Best Answer
In method of characteristics, we reduce the quasilinear PDE as an ODE along characteristic curves and hence solve it for points on a characteristic curve. But at the end of the day, we need to go back and bundle these curves together to form our solution (at least locally).
What you have therefore shown is that, along any characteristic curve, the value of $x+y$ and $xyu$ are constant. But they are allow to change as we move to another characteristic curve (indeed, each characteristic curve is given by $\{x+y=c\}\cap\{xyu=k\}$ locally (some $(c,k)\in\mathbb{R}^2$)). Our solution surface $(x,y,u)$ is locally a 1-parameter family of characteristic curves $(c,k)$ but we don't actually know how $c,k$ are related to each other unless we are given initial data. In other words, we have an arbitrary freedom of reasonable (see next paragraph) parametrisation of this family $(c,k)$ as $\Phi(c,k)=0$, then we have our general solution $\Phi(x+y,xyu)=0$.
Indeed you can check $\Phi(x+y,xyu(x,y))=0$ solves the given PDE locally: differentiating with respect to $x,y$ gives \begin{align*} 0&=\frac{\partial}{\partial x}\Phi(x+y,xyu)\\ &=\Phi_{,1}(x+y,xyu)+\Phi_{,2}(x+y,xyu)(yu+xyu_x)\\ 0&=\frac{\partial}{\partial y}\Phi(x+y,xyu)\\ &=\Phi_{,1}(x+y,xyu)+\Phi_{,2}(x+y,xyu)(xu+xyu_y) \end{align*} If $\Phi_{,2}(x+y,xyu)\neq 0$ (which you must have to say anything about $u$ since the only way for $u$ to come in is through the second slot of $\Phi$), this gives the PDE $$ yu+xyu_x = xu+xyu_y $$ which is just a simple rearrangement of the given PDE. In the general case, we want $\Phi$ to be regular parametrisation (i.e., not both $\Phi_{,1}$ and $\Phi_{,2}$ vanish simultaneously when $\Phi=0$). Similarly for higher dimensions.