An alternative approach to using a matrix representation for the rotation operator is to use clifford algebra instead. You might know clifford algebra in the guise of Pauli and Dirac matrices, from quantum mechanics.
A boost in the $i$th direction is described fully by a "rotor", which takes the form
$$q = \exp(-\gamma_0 \gamma_i \phi/2) = \cosh \frac{\phi}{2} - \gamma_0 \gamma_i \sinh \frac{\phi}{2}$$
If you're used to thinking of the Dirac matrices as, well, matrices, then feel free to write $I \cosh \frac{\phi}{2}$ instead. Here, however, I eschew treating these objects as matrices--I only need their multiplication law to use their clifford algebra properties--and as such, I consider such a term to be "scalar" compared to the vector space formed by the $\gamma_\mu$.
Now, if you're familiar with quaternions, you might recognize what we're doing: indeed, we can proceed exactly as in quaternions (and a good exercise is to derive quaternions using the Pauli algebra). A vector $v$ that is a linear combination $v = v^\mu \gamma_\mu$ can be boosted by
$$R(v) = q v q^{-1} = \exp(-\gamma_0 \gamma_i \phi/2) v \exp(\gamma_0 \gamma_i \phi/2)$$
Multiply all this out, and you get the usual formula for a boost.
Now let's take your example, and let's multiply two rotors: one boosting along $e_x = \gamma_1$ and another boosting along $e_y = \gamma_2$:
$$\left(\cosh \frac{\phi}{2} - \gamma_0 \gamma_1 \sinh \frac{\phi}{2} \right) \left(\cosh \frac{\phi}{2} - \gamma_0 \gamma_2 \sinh \frac{\phi}{2} \right)$$
Remember the clifford multiplication rules: $\gamma_\mu \gamma_\nu = -\gamma_\nu \gamma_\mu$ if $\mu \neq \nu$. and $\gamma_\mu \gamma_\mu = \eta_{\mu \mu}$ (no sum). The gammas are associative, so we can manipulate the expression to get
$$\cosh^2 \frac{\phi}{2} - \left(\sinh\frac{\phi}{2} \cosh \frac{\phi}{2}\right) [\gamma_0 \gamma_1 + \gamma_0 \gamma_2] - \gamma_1 \gamma_2 \sinh^2 \frac{\phi}{2}$$
using $\eta = (+,-,-,-)$ convention.
The presence of the $\gamma_1 \gamma_2$ term tells us that there is a rotation introduced when we compose boosts in this way (at least, when those boosts don't use the same plane).
Does it make sense to talk about the boost velocity for an operation that has both boost terms and pure rotation terms? Remember, this is not merely a boost and then a rotation, nor a rotation and then a boost--they're both happening together as you apply this operator.
If you have a definition of boost velocity that applies in this case, though, I'd be happy to turn the crank and compute it.
As you probably know, there are different irreducible representations of the different symmetry groups that one gets in relativistic quantum field theory. Here we are dealing with the Lorentz symmetry (part of the Poincare symmetry).
One representation, the spin-half representation, is associated with fermions as one finds in the Dirac equation. The indices of this representations are the spin indices, as denoted by the rows and columns of the Dirac matrices. In this representation the Lorentz transformations are represented by spin transformations.
Another representation, the spin-one representation, is associated with vector fields, such as the gauge bosons. It is also the representation for the coordinate vectors. The indices in this case are the space-time indices that are usually denoted by some Greek letter such as $\mu$. In this representation the Lorentz transformations transform the space-time indices (rotations and boosts).
One can combine a pair of spin-half quantities in a why that they act like spin-one quantities. (This is part of a more general process whereby tensor products of irreducible representations of a symmetry group become direct sums of different irreducible representations of that symmetry group.) This means that a product of two spin-half quantities, each transforming by spin transformations, can act like a single spin-one quantity, transforming as a space-time vector, provided that the two quantities are combined in a particular way. It turns out that the Dirac matrices provide the correct way to combine two spin-half quantities so that they transform as a single spin-one quantity. Equation (3.29) demonstrates this property.
Hope this addresses all your questions.
Best Answer
I think the clearest way to think about this is to say that the gamma matrices don't transform. In other words, the fact that they carry a vector index doesn't mean that they form a four vector. This is analogous to how the Pauli matrices work in regular quantum mechanics, so let me talk a little bit about that.
Suppose you have a spin $1/2$ particle in some state $|\psi\rangle$. You can calculate the mean value of $\sigma_x$ by doing $\langle \psi | \sigma_x | \psi\rangle$. Now let's say you rotate your particle by an angle $\theta$ around the $z$-axis. (Warning: There is about a 50% chance my signs are incorrect.) You now describe your particle with a different ket, given by $|\psi'\rangle = \exp(-i \sigma_z \theta /2)$. Remember that we are leaving the coordinates fixed and rotating the system, as is usually done in quantum mechanics. Now the expectation value is given by
$$\langle \psi' | \sigma_x | \psi' \rangle = \langle \psi |\, e^{i\sigma_z \theta /2}\, \sigma_x\, e^{-i \sigma_z \theta / 2}\, | \psi\rangle$$
There is a neat theorem, not too hard to prove, that says that
$$e^{i\sigma_z \theta /2}\, \sigma_x\, e^{-i \sigma_z \theta / 2} = \cos \theta\, \sigma_x -\sin \theta\, \sigma_y$$
So it turns out that the expectation value for the rotated system is also given by $\langle \psi |\, \cos \theta\, \sigma_x -\sin \theta\, \sigma_y \, |\psi\rangle = \cos \theta\, \langle \sigma_x \rangle - \sin \theta\, \langle \sigma_y \rangle$. It's as if we left our particle alone and rotated the Pauli matrices. But note that if we apply the rotation to $|\psi\rangle$, then we don't touch the matrices. Also, I never said that I transformed the matrices. I just transformed the state, and then found out that I could leave it alone and rotate the matrices.
The situation for a Dirac spinor is similar. The analogous identity is that $S(\Lambda) \gamma^\mu S^{-1}(\Lambda) \Lambda^\nu_{\ \mu} = \gamma^\nu$. This is just something that is true; nobody said anything about transforming $\gamma^\mu$. There's no $\gamma^\mu \to \dots$ here.
Now let's take the Dirac equation, $(i \gamma^\mu \partial_\mu - m)\psi = 0$, and apply a Lorentz transformation. This time I will change coordinates instead of boosting the system, but there's no real difference. Let's say we have new coordinates given by $x'^\mu = \Lambda^\mu_{\ \nu} x^\nu$, and we want to see if the Dirac equation looks the same in those coordinates. The field $\psi'$ as seen in the $x'^\mu$ frame is given by $\psi' = S(\Lambda) \psi \iff \psi = S^{-1}(\Lambda) \psi'$, and the derivatives are related by $\partial_\mu = \Lambda^\nu_{\ \mu} \partial'_\nu$. Plugging in we get $(i\gamma^\mu \Lambda^\nu_{\ \mu} \partial'_\nu-m) S^{-1}(\Lambda)\psi' = 0$, which doesn't really look like our original equation. But let's multiply on the left by $S(\Lambda)$. $m$ is a scalar so $S$ goes right through it and cancels with $S^{-1}$. And in the first term we get $S(\Lambda)\gamma^\mu S^{-1}(\Lambda) \Lambda^\nu_{\ \mu}$, which according to our trusty identity is just $\gamma^\nu$. Our equation then simplifies to
$$(i\gamma^\mu \partial'_\mu - m)\psi'=0$$
This is the same equation, but written in the primed frame. Notice how the gamma matrices are the same as before; when you're in class and the teacher writes them on the board, you don't need to ask in what coordinate system they are valid. Everyone uses the same gamma matrices. They're not really a four-vector, but their "transformation law" guarantees that anything written as if they were a four vector is Lorentz invariant as long as the appropiate spinors are present.