This notebook can be used to visualize the IVP for a conservation law PDE:
a(z) z_x + z_y = 0, 
z(x, 0) = f(x)
In applications y is time (t) and we are usually interested in y > 0 (t > 0).

f[x_] := -Tanh[x]
a[z_] := z

pInitialData = ParametricPlot3D[{x, 0, f[x]}, {x, -2, 2}, PlotStyle -> Tube[0.02]]

pSolution = ParametricPlot3D[{x0 + a[f[x0]] y, y, f[x0]}   , {x0, -2, 2}, {y,0,2}, Mesh ->{10, 0}, PlotStyle -> Opacity[.8]]

Show[pInitialData,pSolution, pInitialData, PlotRange -> All, AxesLabel ->{x,y,z}]

characteristics = ParametricPlot[{x0 + a[f[x0]] y, y}   , {x0, -2, 2}, {y,0,2}, Mesh ->{10, 0}, FrameLabel ->{x,y}]

envelope = ParametricPlot[{x0 -((a[f[x0]])/(f'[x0]*a'[f[x0]])), ((-1)/(f'[x0]*a'[f[x0]]))}, {x0, -1,1}, PlotStyle -> Red]

Show[characteristics,envelope]

Now I will do the same thing for a traffic flow problem.  I have also shown only the final figures of each set.  The parts of the figure are defined but not plotted with a “;” and then the composite figure is shown.
Here the characteristics are NOT the paths of individual cars, but the lines on which the density is constant.
It takes some getting used to.
All the cars are moving to the right, and normalized density of cars is z \in [0, 1].

f[x_] := (2 + Tanh[x])/3
a[z_] := 1-2z
pInitialData = ParametricPlot3D[{x, 0, f[x]}, {x, -2, 2}, PlotStyle -> Tube[0.02]];
pSolution = ParametricPlot3D[{x0 + a[f[x0]] y, y, f[x0]}   , {x0, -2, 2}, {y,0,3}, Mesh ->{10, 0}, PlotStyle -> Opacity[.8]];
Show[pInitialData,pSolution, pInitialData, PlotRange -> All, AxesLabel ->{x,y,z}]
characteristics = ParametricPlot[{x0 + a[f[x0]] y, y}   , {x0, -2, 2}, {y,0,3}, Mesh ->{10, 0}, FrameLabel ->{x,y}];
envelope = ParametricPlot[{x0 -((a[f[x0]])/(f'[x0]*a'[f[x0]])), ((-1)/(f'[x0]*a'[f[x0]]))}, {x0, -1,1}, PlotStyle -> Red];
(* To check the position of the shock (computed by hand), plot the point below *)
Show[characteristics,envelope, Graphics@Point[{-1/2, 3/2}]]