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.
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.
Best Answer
OK I think I see where your confusion lies.
You're talking about the four velocity and acceleration in the instantaneous rest frame $F'$, and as you say in this frame the four velocity is $(1, 0, 0, 0)$. Your mistake is to assume the four velocity is constant in $F'$, because it is not. Remember that after an infinitesimal time $dt$ the rocket is not longer in $F'$ - it is in a new instantaneous rest frame $F''$. The rocket's velocity in the new rest frame $F''$ is still $(1, 0, 0, 0)$, but in old $F'$ frame it has now changed due to the acceleration. Hence $d{\bf u}/dt$ in $F'$ is not zero.
If you're interested chapter 6 of Gravitation by Misner, Thorne and Wheeler derives the equations of motion that you started with.