We know more than just the fact that $M$ has eigenvalues in $\{-2,-1,0,1\}$. If $M$'s eigenvalues are $\lambda_1, \lambda_2$, we know that $M^2+M$ has eigenvalues $\lambda_1^2+\lambda_1$ and $\lambda_2^2+\lambda_2$. Since $M^2+M$ has eigenvalues $0$ and $2$, we know that one of $M$'s eigenvalues (say, $\lambda_1$) has $\lambda_1^2+\lambda_1 = 0$ and the other eigenvalue has $\lambda_2^2+\lambda_2 = 2$.
So one eigenvalue is $0$ or $-1$, and the other eigenvalue is $-2$ or $1$.
All four combinations are possible. Write the matrices in the basis $\mathcal B$ of the eigenvectors $(1,1)$ and $(1,-1)$: in this basis, we have
$$
(M_{\mathcal B})^2 + M_{\mathcal B} = \begin{bmatrix}2 & 0 \\ 0 & 0\end{bmatrix}
$$
and so any of the possibilities
$$
M_{\mathcal B} \in \left\{\begin{bmatrix}-2 & 0 \\ 0 & 0\end{bmatrix}, \begin{bmatrix}-2 & 0 \\ 0 & -1\end{bmatrix}, \begin{bmatrix}1 & 0 \\ 0 & 0\end{bmatrix}, \begin{bmatrix}1 & 0 \\ 0 & -1\end{bmatrix}\right\}
$$
will work. Convert these to the standard basis and you find four possibilities for $M$, each corresponding to one of the four cases for $M$'s two eigenvalues.
(It's worth spelling out in detail why this argument still makes sense without assuming that $M$ is diagonalizable. We don't need to assume $\lambda_1 \ne \lambda_2$ for the statement about the eigenvalues of $M^2+M$ to hold: if $\lambda_1 = \lambda_2$ and this eigenvalue has algebraic multiplicity $2$, then $M^2+M$ has $\lambda_1^2+\lambda_1 = \lambda_2^2+\lambda_2$ as an eigenvalue with algebraic multiplicity $2$. This is then contradicted by our knowledge of the eigenvalues of $M^2+M$.)
Set the components of the matrix equal to one another.
On one hand, the product of the two boosts is equal to:
$$\begin{bmatrix}\gamma_1 \gamma_2 (1+\beta_1\beta_2) & - \gamma_1\gamma_2 (\beta_1 + \beta_2) \\ - \gamma_1\gamma_2 (\beta_1 + \beta_2) & \gamma_1 \gamma_2 (1+\beta_1\beta_2)\end{bmatrix}$$
on the other hand, for closure, you want this product to have the form of a boost for some unknown value of $\beta_3$ and $\gamma_3(\beta_3)$:
$$\begin{bmatrix}\gamma_3 & -\gamma_3 \beta_3 \\ -\gamma_3\beta_3 & \gamma_3\end{bmatrix}$$
Setting these two matrices equal to one another, you find constraints on $\gamma_3$ and $\beta_3$:
$$\gamma_3 = \gamma_1\gamma_2(1+\beta_1\beta_2)\\ -\gamma_3\beta_3 = -\gamma_1\gamma_2(\beta_1+\beta_2) \\ \Rightarrow \beta_3 = \frac{\beta_1+\beta_2}{1+\beta_1\beta_2}$$
So you can solve for $\gamma_3$ and $\beta_3$ in terms of the original parameters $\beta_1,\beta_2,\gamma_1,\gamma_2$ and show that $\gamma_3$ is indeed the Lorenz factor corresponding to $\beta_3$, i.e. you can prove that, as required, $$\gamma_3 = \frac{1}{\sqrt{1-\beta_3^2}}$$
so the product of two boosts is itself a boost corresponding to $\beta_3,\gamma_3$.
Best Answer
The matrix $$N=\begin{pmatrix} 0 & 1\\ 0 & 0 \end{pmatrix}$$ is nilpotent with index 2 of nilpotency: $N^2=0$ so by the binomial formula we have
$$\begin{pmatrix} 1 & 1\\ 0 & 1 \end{pmatrix}^n=(I_2+N)^n=\sum_{k=0}^n {n\choose k}N^k={n\choose 0}I_2+{n\choose 1}N=I_2+nN=\begin{pmatrix} 1& n\\ 0 & 1 \end{pmatrix}$$