Solution (source code)

= Solution

The <method of characteristics> uses
$$
\frac{dx}{ds}=a(x,y,u),\qquad\frac{dy}{ds}=b(x,y,u),\qquad\frac{du}{ds}=c(x,y,u).
$$
For a graph parametrized by $(x,y,u(x,y))$, a <normal vector> is $N=(-u_x,-u_y,1)$. The PDE says $N\cdot(a,b,c)=0$, so the characteristic vector is tangent to the solution surface. Flowing a suitable initial curve along these directions constructs a surface; where its projection to the $(x,y)$ plane has nonzero <Jacobian determinant>, the <inverse function theorem> makes it a graph. Local uniqueness of this general flow needs the usual regularity assumptions beyond mere continuity. The particular characteristic system below is smooth and can be solved explicitly.

Label the initial points by $\eta\ge0$. Since $dx/ds=1$, take $x=s$. The other two equations integrate to
$$
u=1+\sin\eta+\frac{x^2}{2},\qquad
y=\eta+x(1+\sin\eta)+\frac{x^3}{6}.
$$
Eliminating $\eta=y-xu+x^3/3$ gives the implicit solution
$$
\boxed{u=1+\frac{x^2}{2}+\sin\!\left(y-xu+\frac{x^3}{3}\right).}
$$
For each fixed $|x|<1$, the <derivative> of the displayed $y$ with respect to $\eta$ is $1+x\cos\eta\ge1-|x|>0$. It increases from its value at $\eta=0$ to infinity, so it is a bijection from $[0,\infty)$ onto $[y_L(x),\infty)$, where
$$
\boxed{y_L(x)=x+\frac{x^3}{6}.}
$$
Thus there is exactly one characteristic label, and hence one solution value, throughout the requested region, including the initial-label boundary as a trace.

The <characteristic fold for sinusoidal data in a forced Burgers equation> is seen by differentiating with respect to $y$:
$$
u_y=\frac{\cos\eta}{1+x\cos\eta}.
$$
For a fixed <integer> $k\ge0$, let $y_*=(2k+1)\pi+7/6$. As $x\to1^-$ at this fixed $y_*$, the label tends to $(2k+1)\pi$, giving $u_y\to-\infty$. To justify the label limit, the limiting map at $x=1$ is strictly increasing despite isolated zero <derivatives>, so its inverse has a unique continuous limit. Similarly, at any fixed $y_*=2k\pi-7/6$ with $k\ge1$, the limit as $x\to-1^+$ has $\eta\to2k\pi$ and $u_y\to+\infty$. These fixed values lie within the region near the respective boundary, proving the requested gradient blow-up.