In my previous answer, I mentioned that the form of a general boost to speed $\beta c$ in the direction of the unit vector $\hat n$ is
$$B(\beta,\hat n)=
\left(
\begin{array}{cccc}
\gamma & -\gamma \beta n_x & -\gamma \beta n_y & -\gamma \beta n_z \\
-\gamma \beta n_x & 1+(\gamma-1)n_x^2 & (\gamma-1)n_xn_y & (\gamma-1)n_xn_z \\
-\gamma \beta n_y & (\gamma-1)n_yn_x & 1+(\gamma-1)n_y^2 & (\gamma-1)n_yn_z \\
-\gamma \beta n_z & (\gamma-1)n_zn_x & (\gamma-1)n_zn_y & 1+(\gamma-1)n_z^2 \\
\end{array}
\right)
$$
where
$$\gamma=\frac{1}{\sqrt{1-\beta^2}}.$$
A general rotation by angle $\theta$ around an axis determined by the unit vector $\hat u$ (with rotational direction given by the right-hand-rule) is
$R(\theta,\hat u)=$
$$
\left(
\begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 &
\cos \theta +u_x^2 \left(1-\cos \theta\right) &
u_x u_y \left(1-\cos \theta\right) - u_z \sin \theta &
u_x u_z \left(1-\cos \theta\right) + u_y \sin \theta \\
0 &
u_y u_x \left(1-\cos \theta\right) + u_z \sin \theta &
\cos \theta + u_y^2\left(1-\cos \theta\right) &
u_y u_z \left(1-\cos \theta\right) - u_x \sin \theta \\
0 &
u_z u_x \left(1-\cos \theta\right) - u_y \sin \theta &
u_z u_y \left(1-\cos \theta\right) + u_x \sin \theta &
\cos \theta + u_z^2\left(1-\cos \theta\right)
\end{array}
\right)
$$
By multiplying these matrices, you can see that a general boost followed by a general rotation has the form
$$R(\theta,\hat u)B(\beta,\hat n)=
\left(
\begin{array}{cccc}
\gamma & -\gamma\beta n_x & -\gamma\beta n_y & -\gamma\beta n_z \\
- & - & - & - \\
- & - & - & - \\
- & - & - & - \\
\end{array}
\right)
$$
and a general rotation followed by a general boost has the form
$$B(\beta,\hat n)R(\theta,\hat u)=
\left(
\begin{array}{cccc}
\gamma & - & - & - \\
-\gamma\beta n_x & - & - & - \\
-\gamma\beta n_y & - & - & - \\
-\gamma\beta n_z & - & - & - \\
\end{array}
\right)
$$
Here the matrix elements indicated by a dash are complicated expressions involving both the boost parameters and the rotation parameters. But the first row (in the $RB$ case) or the first column (in the $BR$ case) is simple: It depends only on the boost parameters!
This means that if we want to break down a general Lorentz transformation matrix into a boost and a rotation, we can just extract the boost parameters from the first row or column.
Let's see how this works in the example I gave of composing two boosts, $B\left(\frac12, \hat x\right)$ followed by $B\left(\frac12, \hat y\right)$. Multiplying the two boost matrices gives
$$\begin{align}
B\left(\frac12, \hat y\right)B\left(\frac12, \hat x\right)&=
\left(
\begin{array}{cccc}
\frac{2}{\sqrt{3}} & 0 & -\frac{1}{\sqrt{3}} & 0 \\
0 & 1 & 0 & 0 \\
-\frac{1}{\sqrt{3}} & 0 & \frac{2}{\sqrt{3}} & 0 \\
0 & 0 & 0 & 1 \\
\end{array}
\right)
\left(
\begin{array}{cccc}
\frac{2}{\sqrt{3}} & -\frac{1}{\sqrt{3}} & 0 & 0 \\
-\frac{1}{\sqrt{3}} & \frac{2}{\sqrt{3}} & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
\end{array}
\right)\\
&=
\left(
\begin{array}{cccc}
\frac{4}{3} & -\frac{2}{3} & -\frac{1}{\sqrt{3}} & 0 \\
-\frac{1}{\sqrt{3}} & \frac{2}{\sqrt{3}} & 0 & 0 \\
-\frac{2}{3} & \frac{1}{3} & \frac{2}{\sqrt{3}} & 0 \\
0 & 0 & 0 & 1 \\
\end{array}
\right)
\end{align}
$$
The matrix on the right (the composition of the two boosts) is clearly not a boost, because it isn't symmetric. And it's clearly not a rotation, because it mixes space and time coordinates. It is a more general Lorentz transformation that is neither a boost nor a rotation but a combination of them.
If we want to write it in $RB$ (boost-then-rotation) form, we extract the boost parameters from the first row. From the top-left element we find that
$$\gamma=\frac43$$
and thus the boost speed is
$$\beta=\sqrt{1-\frac{1}{\gamma^2}}=\frac{\sqrt7}{4}$$
and the product $\gamma\beta$ is
$$\gamma\beta=\frac{\sqrt{7}}{3}.$$
From the other three elements in the top row, we can easily find the boost direction $\hat n$. We have
$$-\gamma\beta n_x=-\frac23$$
$$-\gamma\beta n_y=-\frac{1}{\sqrt{3}}$$
$$-\gamma\beta n_z=0$$
and thus
$$\hat n=\left(\frac{2}{\sqrt{7}},\frac{\sqrt{3}}{\sqrt{7}},0\right).$$
Once you have extracted the parameters $\beta$ and $\hat n$ for the boost part of a general Lorentz transformation $\Lambda=RB$ from the first row of $\Lambda$, you can use the general boost formula to compute $B$. In our case, it is
$$\begin{align}
B&=B\left(\frac{\sqrt{7}}{4},\left(\frac{2}{\sqrt{7}},\frac{\sqrt{3}}{\sqrt{7}},0\right)\right)\\
&=
\left(
\begin{array}{cccc}
\frac{4}{3} & -\frac{2}{3} & -\frac{1}{\sqrt{3}} & 0 \\
-\frac{2}{3} & \frac{25}{21} & \frac{2}{7 \sqrt{3}} & 0 \\
-\frac{1}{\sqrt{3}} & \frac{2}{7 \sqrt{3}} & \frac{8}{7} & 0 \\
0 & 0 & 0 & 1 \\
\end{array}
\right)
\end{align}
$$
One can then invert it and find the rotation factor $R=\Lambda B^{-1}$.
If you prefer to write a general transformation in the form $\Lambda=BR$, then extract the boot parameters from the first column of $\Lambda$, compute $B$, and then compute $R=B^{-1}\Lambda$.
In either case, computing the inverse boost factor $B^{-1}$ is easy: Rather than actually inverting a $4\times 4$ matrix, you can just reverse the direction of the boost parameter $\hat n$. So in our example,
$$\begin{align}
B^{-1}&=\left[B\left(\frac{\sqrt{7}}{4},\left(\frac{2}{\sqrt{7}},\frac{\sqrt{3}}{\sqrt{7}},0\right)\right)\right]^{-1}\\
&=B\left(\frac{\sqrt{7}}{4},\left(-\frac{2}{\sqrt{7}},-\frac{\sqrt{3}}{\sqrt{7}},0\right)\right)\\
&=
\left(
\begin{array}{cccc}
\frac{4}{3} & \frac{2}{3} & \frac{1}{\sqrt{3}} & 0 \\
\frac{2}{3} & \frac{25}{21} & \frac{2}{7 \sqrt{3}} & 0 \\
\frac{1}{\sqrt{3}} & \frac{2}{7 \sqrt{3}} & \frac{8}{7} & 0 \\
0 & 0 & 0 & 1 \\
\end{array}
\right)
\end{align}$$
The rotation matrix $R$ satifying $\Lambda=RB$ is then
$$\begin{align}
R&=\Lambda B^{-1}\\
&=
\left(
\begin{array}{cccc}
\frac{4}{3} & -\frac{2}{3} & -\frac{1}{\sqrt{3}} & 0 \\
-\frac{1}{\sqrt{3}} & \frac{2}{\sqrt{3}} & 0 & 0 \\
-\frac{2}{3} & \frac{1}{3} & \frac{2}{\sqrt{3}} & 0 \\
0 & 0 & 0 & 1 \\
\end{array}
\right)
\left(
\begin{array}{cccc}
\frac{4}{3} & \frac{2}{3} & \frac{1}{\sqrt{3}} & 0 \\
\frac{2}{3} & \frac{25}{21} & \frac{2}{7 \sqrt{3}} & 0 \\
\frac{1}{\sqrt{3}} & \frac{2}{7 \sqrt{3}} & \frac{8}{7} & 0 \\
0 & 0 & 0 & 1 \\
\end{array}
\right)\\
&=
\left(
\begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 & \frac{4 \sqrt{3}}{7} & -\frac{1}{7} & 0 \\
0 & \frac{1}{7} & \frac{4 \sqrt{3}}{7} & 0 \\
0 & 0 & 0 & 1 \\
\end{array}
\right)
\end{align}
$$
This is recognizable as a rotation around the $z$-axis because it mixes only $x$ and $y$ together. Rotations by $\theta$ around the $z$-axis look like
$$
\left(
\begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 & \cos\theta & -\sin\theta & 0 \\
0 & \sin\theta & \cos\theta & 0 \\
0 & 0 & 0 & 1 \\
\end{array}
\right)
$$
so we immediately see that the rotation angle is such that
$$\cos\theta=\frac{4\sqrt{3}}{7}$$
and
$$\sin\theta=\frac{1}{7}.$$
But this raises the question of how to extract the rotation angle and rotation direction if $R$ doesn't have such a simple form as in this example.
By taking the trace of the general rotation matrix, we see that
$$\text{Tr}R(\theta,\hat u)=2+2\cos\theta$$
so the rotation angle of a general rotation $R$ is
$$\theta=\cos^{-1}\frac{\text{Tr}R-2}{2}.$$
And, since the unit vector $\hat u$ along the rotation axis doesn't get rotated by the rotation,
$$R\left(\begin{array}{c}1\\u_x\\u_y\\u_z\end{array}\right)
=\left(\begin{array}{c}1\\u_x\\u_y\\u_z\end{array}\right)$$
This says that $\hat u$ is the spatial part of an eigenvector of $R$ with eigenvalue $1$.
Let's see how this works in an example where $R$ is a more complicated rotation. Consider the Lorentz transformation
$$\Lambda=
\left(
\begin{array}{cccc}
\frac{8}{3 \sqrt{3}} & -\frac{4}{3 \sqrt{3}} & -\frac{2}{3} & -\frac{1}{\sqrt{3}} \\
-\frac{1}{\sqrt{3}} & \frac{2}{\sqrt{3}} & 0 & 0 \\
-\frac{2}{3} & \frac{1}{3} & \frac{2}{\sqrt{3}} & 0 \\
-\frac{4}{3 \sqrt{3}} & \frac{2}{3 \sqrt{3}} & \frac{1}{3} & \frac{2}{\sqrt{3}} \\
\end{array}
\right)
$$
which I constructed by composing three boosts: the first by $c/2$ along $\hat x$, the second by $c/2$ along $\hat y$, and the third by $c/2$ along $\hat z$.
Extracting the boost parameters as before from the first row, we find
$$\beta=\frac{\sqrt{37}}{8}$$
and
$$\hat n=\frac{(4,2\sqrt{3},3)}{\sqrt{37}}.$$
After working out $B$ and $B^{-1}$, one finds the rotation to be
$$R=
\left(
\begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 & \frac{2}{37} \left(6+7 \sqrt{3}\right) & \frac{2}{37} \left(3 \sqrt{3}-8\right) & \frac{1}{37} \left(9-8 \sqrt{3}\right) \\
0 & \frac{1}{37} \left(8 \sqrt{3}-9\right) & \frac{2}{37} \left(6+7 \sqrt{3}\right) & \frac{2}{37} \left(3 \sqrt{3}-8\right) \\
0 & -\frac{2}{37} \left(3 \sqrt{3}-8\right) & \frac{1}{37} \left(8 \sqrt{3}-9\right) & \frac{2}{37} \left(6+7 \sqrt{3}\right) \\
\end{array}
\right)
$$
The trace formula gives the rotation angle as
$$\theta=\cos^{-1}\frac{42\sqrt{3}-1}{74}\approx 14.18^\circ.$$
An unnormalized spatial eigenvector with eigenvalue $1$ is $(0,1,-1,1)$. (I suggest using a computer algebra system to compute the eigenvalues and eigenvectors! There is also a trivial temporal eigenvector $(1,0,0,0)$ also with eigenvalue $1$, and two complex eigenvectors with complex eigenvalues $e^{i\theta}$ and $e^{-i\theta}$. The sum of the eigenvalues equals the trace value, $2+2\cos\theta$, as it should.)
Therefore the rotation axis is the normalized spatial part of this eigenvector,
$$\hat u=\frac{(1,-1,1)}{\sqrt{3}}.$$
