starting with 2 particle collision :
\begin{align*}
&m_1\,(\mathbf v_1-\mathbf u_1)=-\lambda_{12}\,\mathbf{\hat{r}}_{12}\\
&m_2\,(\mathbf v_2-\mathbf u_2)=+\lambda_{12}\,\mathbf{\hat{r}}_{12}\\
&\left[~\mathbf v_2-\mathbf v_1+\epsilon(\mathbf u_2-\mathbf u_1)~\right]\cdot \,\mathbf{\hat{r}}_{12}=0\\
&\text{with}\\
&\mathbf{\hat{r}}_{12}=\frac{\mathbf r_2-\mathbf r_1}{\parallel\mathbf r_2-\mathbf r_1\parallel}
\end{align*}
where $~\mathbf v_i~$ the velocity after the collision, $~\mathbf u_i~$ the starting velocity and $~\epsilon~$ the restoration coefficient.
$~(\epsilon=1~$ elastic collision).
you obtain 5 equations for the 5 unknowns $~\mathbf v_1~,\mathbf v_2~,\lambda_{12}~$
writing the equations with matrix notation $~\mathbf A\,\mathbf x=\mathbf b$
\begin{align*}
&\underbrace{\begin{bmatrix}
m_1\,\mathbf I_2 & 0\mathbf I_2 &+\mathbf{\hat{r}}_{12}^T \\
0\mathbf I_2 & m_2\,\mathbf I_2 & -\mathbf{\hat{r}}_{12}^T \\
-\mathbf{\hat{r}}_{12}^T & +\mathbf{\hat{r}}_{12}^T & 0 \\
\end{bmatrix}}_{\mathbf A_{5\times 5}}
\underbrace{ \begin{bmatrix}
\mathbf v_1 \\
\mathbf v_2 \\
\lambda_{12} \\
\end{bmatrix}}_{\mathbf x}=
\underbrace{\begin{bmatrix}
m_1\,\mathbf u_1 \\
m_2\,\mathbf u_2 \\
\epsilon\,\mathbf{\hat{r}}_{12}\cdot (\mathbf u_1-\mathbf u_2) \\
\end{bmatrix}}_{\mathbf b}\tag 1
\end{align*}
you can adapt equation (1) for collision between any two particle, for example three particle collision, the possible collision are between
particle 1 and 2, 1 and 3 and 2 and 3,for each pair you have to solve the matrix equation (1).