Here is a sketch of the characteristic curves in the $x$-$t$ plane, which equation is given by
$$
\begin{aligned}
x'(t) &= u(x(t),t) \\
&= u(x(0),0)
\end{aligned}
$$
Those curves intersect already at time zero:
Thus, a shock wave arises.
The Burgers' equation is rewritten in the conservative form $u_t + f(u)_x = 0$, where $f(u) = \frac{1}{2}u^2$. The shock wave with speed $s$, left state $u_L = 1$ and right state $u_R = 0$ writes
$$
u(x,t) =
\left\lbrace
\begin{aligned}
&1 &&\text{if}\quad x<st \\
&0 &&\text{if}\quad st<x \, .
\end{aligned}
\right.
$$
The speed of shock must satisfy the Rankine-Hugoniot condition $s = \frac{f(u_R) - f(u_L)}{u_R - u_L}$, i.e.
$s = \frac{1}{2}$.
Here is a modified sketch of the $x$-$t$ plane, which accounts for the shock:
After the basic shock waves $x=t/2$ and $x=t+2 $ (shown in red) are accounted for, the characteristics look like this.
This picture is correct for times $t<1$, but needs two adjustments after that. Note that the solution within rarefaction wave is $u(x,t) = (x-1)/t$.
Rarefaction catches up with shock on the right: $x=3$, $t=1$.
From this point on, the position of this shock satisfies
$$
\frac{dx}{dt} = \frac12 \left( \frac{x-1}{t} + 0 \right) = \frac{x-1}{2t}
$$
hence $x(t) = 2\sqrt{t} +1$.
Shock on the left catches up with rarefaction: $x=1$, $t=2$.
From this point on, the position of this shock satisfies
$$
\frac{dx}{dt} = \frac12 \left( 1 + \frac{x-1}{t} \right) = \frac{x-1+t}{2t}
$$
hence $x(t) = t-\sqrt{2t}+1$.
Notice that in both cases the space-time trajectory of the shock becomes parabolic: this always happens when a rarefaction wave interacts with a constant velocity field.
Plotting the new trajectories of shocks, and terminating the characteristics accordingly, we get the following picture (one has to look closely to notice the curvature of red curves).
Final battle: two shocks come together
At the time $t=6+4\sqrt{2} \approx 11.5$ (off the chart above), the two shock waves meet at $x=5+2\sqrt{2}$. At this moment, rarefaction ceases to exist. The single shock wave forms, separating the velocities $1$ and $0$. It propagates according to the equation
$$
x(t) = 5+2\sqrt{2} + \frac12 (t-6-4\sqrt{2}),\quad t\ge 6+4\sqrt{2}
$$
Best Answer
If you want a numerical way of finding the shock, you could take your array $u[i,j]$ and consider the partial derivative in the x-direction:
dudx = diff(u) / dx
And then if you contour plot dudx, you should get an image where the shock shows up as yellow. This should be along the line $t = 2x$, which you can see in your plot above.