  We prove a stable shock formation result for a large class of systems of
quasilinear wave equations in two spatial dimensions. We give a precise
description of the dynamics all the way up to the singularity. Our main theorem
applies to systems of two wave equations featuring two distinct wave speeds and
various quasilinear and semilinear nonlinearities, while the solutions under
study are (non-symmetric) perturbations of simple outgoing plane symmetric
waves. The two waves are allowed to interact all the way up to the singularity.
Our approach is robust and could be used to prove shock formation results for
other related systems with many unknowns and multiple speeds, in various
solution regimes, and in higher spatial dimensions. However, a fundamental
aspect of our framework is that it applies only to solutions in which the
"fastest wave" forms a shock while the remaining solution variables do not.
  Our approach is based on an extended version of the geometric vectorfield
method developed by D. Christodoulou in his study of shock formation for scalar
wave equations as well as the framework developed in our recent joint work with
J. Luk, in which we proved a shock formation result for a quasilinear
wave-transport system featuring a single wave operator. A key new difficulty
that we encounter is that the geometric vectorfields that we use to commute the
equations are, by necessity, adapted to the wave operator of the
(shock-forming) fast wave and therefore exhibit very poor commutation
properties with the slow wave operator, much worse than their commutation
properties with a transport operator. To overcome this difficulty, we rely on a
first-order reformulation of the slow wave equation, which, though somewhat
limiting in the precision it affords, allows us to avoid uncontrollable
commutator terms.
