How do we solve Burgers' equation using the method of characteristics? I want to solve \begin{align}&\forall(t,x)\in(0,\infty)\times\mathbb R:\left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}\right)(t,x)=0;\tag1\\&\forall x\in\mathbb R:u(0,x)=u_0(x)\tag2,\end{align} where $\Omega:=(0,\infty)\times\mathbb R$ and $u_0\in C^1(\mathbb R)$, using the method of characteristics. Assume that $u_0'\le0$ and $$u_0(x)=\left.\begin{cases}1&\text{, if }x\le0\\0&\text{, if }x\ge1\end{cases}\right\}\;\;\;\text{for all }x\in\mathbb R\tag3.$$
As usual, we consider a path $\gamma$ with $\gamma'=(1,\underbrace{u\circ\gamma}_{=:\:z})$, since then $(1)$ implies $z'=0$ and hence $z=c_1$ for some $c_1\in\mathbb R$, which in turn impleis $\gamma=(c_2+t,c_3+c_1t)$ for some $c_2,c_3\in\mathbb R$.
Fix $(t_0,x_0)\in(0,\infty)\times\mathbb R$. Since we are interested in the value $u(t_0,x_0)$, we choose $c_2:=0$ and $c_3:=x_0-c_1t_0$ to obtain $$u(t_0,x_0)=z(t_0)=z(0)=u_0(x_0-c_1t_0)\tag4.$$

Question 1: So, is this the whole story? Do we obtain different values for the solution of $(1)$ by varying $c_1$? And does this that $(1)$ and $(2)$ together do not admit a unique solution? (And do they admit a solution at all?)
Question 2: By $(3)$, we see that if $c_1\ge\frac{x_0}{t_0}$, then $u(t_0,x_0)=1$ ad if $c_1\le\frac{x_0-1}{t_0}$, then $u(t_0,x_0)=0$, but what can we infer from that? Does this somehow yield a contradiction for general $t_0$? Or, phrased differently, does this imply that we won't be able to obtain a "global" solution, but only a solution up to some finite time $T>0$?

 A: In the sense of strong solutions (functions satisfying the differential equation as written), the PDE is not satisfied where the derivative does not exist. For this problem, the initial condition has no derivative at x=0. The strong solution is well defined everywhere in $(x,t)$ except along the line $x=t$, due to the localization along characteristics of hyperbolic PDEs. In this sense, there are no global strong solutions. However, there are global weak solutions.
Example:
To investigate further, even smooth initial conditions lead to discontinuities forming for $t>0$, at the intersection of characteristics.
Denoting the characteristic variable as $\xi$, one has
$\frac{dx}{dt} = u_0(\xi)\implies\xi = x - u_0(\xi)t$,
and differentiating with respect to $x$ gives $\xi_x = (1+t u_0'(\xi))^{-1}$,
so that the derivative ceases to exist at the critical times $t^*(\xi) = -(u_0'(\xi))^{-1}$.
For example, consider
$u_0(0,x) = -\text{arctan}(x)\implies t^*(\xi) = 1+\xi^2$.
The earliest such time of shock formation occurs for $\xi=0$ at $t=1$.
Weak solutions: Since the PDE does not allow us to determine behavior along characteristics beyond the formation of shocks, one has to consider the weak form of the PDE,
$u_t + uu_x = 0\quad\implies\quad \frac{d}{dt}\int_{x_1}^{x_2}udx + \frac{1}{2}u^2\Big|_{x_1}^{x_2} = 0$
by integrating over a finite region $x\in [x_1, x_2]$. Provided that the derivatives are well defined, then the limit $x_2\to x_1$ recovers the local form of the PDE as originally posed. However, the weak or integral form of the equation is still defined even where derivatives do not exist. This allows us to construct global solutions. If you think about it, most problems in physics are actually integral equations (conservation relations) which are reduced to their point differential equation forms by supposing smooth solutions. So really, we've "gone backwards" by starting at the point-wise PDE.
To see that the weak solution has a well-defined speed at discontinuities, consider a point $x=s(t)$ of discontinuity,
$\frac{d}{dt}\Big(\int_{x_1}^{s(t)}udx + \int_{s(t)}^{x_2}udx\Big) + \frac{1}{2}u^2\Big|_{x_1}^{x_2} = 0$
Using the Leibnitz rule,
$\frac{d}{dt}\int_{a(t)}^{b(t)} u dx = b'(t)u(a_{-}) - a'(t)u(b_{+}) + \int_{a(t)}^{b(t)}u_tdx$
the weak form is equivalently
$s'(t)(u(x_1) - u(x_2)) + \int_{x_1}^{x_2} u_t dx = \frac{1}{2}(u^2(x_1) - u^2(x_2))$
so in the limit as $x_1$ and $x_2$ approach $s(t)$ from their respective sides, the speed of the point $s(t)$ is
$s'(t) = \frac{1}{2}\frac{u_{-}^2 - u_{+}^2}{u_{-} - u_{+}} = \frac{1}{2}(u_{-} + u_{+})$
or the average of the two sides of the discontinuity. For more general conservation laws this is called the Rankine-Hugoniot condition: the shock speed is the jump in the flux divided by the jump in the conserved quantity.
This shows that the weak form has a solution even at points of discontinuity, leading to global solutions of the problem.
A: Let's outline a general method for first order quasilinear PDEs in two variables coupled with an initial condition specified along some curve:
$$\begin{cases}
a(x,y,u(x,y))~\partial_xu(x,y)+b(x,y,u(x,y))~\partial_yu(x,y)=c(x,y,u(x,y))\\
u(\gamma_1(r),\gamma_2(r))=\phi(r)\\
\end{cases}$$
We make a slight generalization to the normal characteristic equations and consider the coupled PDEs
$$ \begin{array}{l}
\partial _{s} x( r,s) =a( x( r,s) ,y( r,s) ,z( r,s))\\
\partial _{s} y( r,s) =b( x( r,s) ,y( r,s) ,z( r,s))\\
\partial _{s} z( r,s) =c( x( r,s) ,y( r,s) ,z( r,s))
\end{array}$$
Subject to initial conditions
$$ \begin{array}{l}
x( r,0) =\gamma _{1}( r)\\
y( r,0) =\gamma _{2}( r)\\
z( r,0) =\phi ( r)
\end{array}$$
Once we have found solutions for $x(r,s)~,~y(r,s)~,~z(r,s)$ we first invert these three equations to find expressions for $r,s$ in terms of $x,y$ and then use $$z(r,s)=u(x(r,s),y(r,s))$$
Hopefully the notation isn't too confusing. I think things will be simpler when we have a concrete example. We have the inviscid Burgers' equation coupled with a linear initial condition:
$$\begin{cases}
u( x,t) \ \partial _{x} u( x,t)+\partial _{t} u( x,t)  =0\\
u( r,0) =\phi(r)=k_{1} r+k_{2}
\end{cases}$$
So our curve along which we specify the initial conditions is
$$\gamma(r)=(\gamma_1(r),\gamma_2(r))=(r,0)$$
And our functions $a,b,c$ are
$$a(x_1,x_2,x_3)=x_3~~,~~b=1~~,~~c=0$$
So, our coupled PDEs are
$$ \begin{array}{l}
\partial _{s} x( r,s) =z( r,s)\\
\partial _{s} t( r,s) =1\\
\partial _{s} v( r,s) =0
\end{array}$$
With initial conditions
$$ \begin{array}{l}
x( r,0) =r\\
t( r,0) =0\\
v( r,0) =k_{1} r+k_{2}
\end{array}$$
Which leads us to solutions (work this out for yourself)
$$ \begin{array}{l}
x( r,s) =k_{1} rs+k_{2} s+r\\
t( r,s) =s\\
v( r,s) =k_{1} r+k_{2}
\end{array}$$
We invert the first two equations to yield
$$r(x,t)=\frac{x-k_2t}{1+k_1t}~~,~~s(x,t)=t$$
So therefore
$$v(r(x,t),s(x,t))=u(x,t)=k_1~\frac{x-k_2t}{1+k_1t}+k_2=\frac{k_1x+k_2}{k_1t+1}$$
One can check (preferably with a CAS) that this satisfies our PDE and the prescribed initial conditions. In fact, this result is already known due to Chandrasekhar in 1943. Coming back to your question, you had an initial condition
$$u(r,0)=1-\operatorname{step}(r)$$
Where $\operatorname{step}$ indicates the Heaviside step function, which unfortunately is not linear, but, we can split it into linear pieces. On the positive real axis this means we have
$$u_{\text{right}}(x,t)=0$$
Similarly
$$u_{\text{left}}(x,t)=1$$
So, if we are unwilling to convert the PDE into its weak form as proposed by @bartholovidus then this is the only possible solution we can have. But, this solution is undefined for $x=0$, since the initial condition is discontinuous there. In this sense no global solutions exist.
