Here goes my own answer (I will leave it unaccepted until someone comes up with a better, more elegant, more efficient, solution). As per @Jean-ClaudeArbaut's suggestion in the comments (which however overlooked the $m_i$ should be naturals), I express the matrix in row echelon form and choose 4 parameters (the difference between the dimension [6] and rank [2] of my matrix):
$$
\begin{pmatrix}
n_1 \\
n_2 - 2n_1 \\
m_3 \\
m_4 \\
m_5 \\
m_6
\end{pmatrix}
=
\begin{pmatrix}
1 & 1 & 0 & 0 & 1 & 0 \\
0 & -1 & 2 & 3 & -2 & 1 \\
0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
m_1 \\
m_2 \\
m_3 \\
m_4 \\
m_5 \\
m_6
\end{pmatrix}
$$
I have chosen $m_3$ to $m_6$ as parameters out of convenience because the two rows of my original matrix allow me to easily establish the bounds to find solutions:
$$
m_3 \in [0,\text{int}(n_2/2)] \\
m_4 \in [0,\text{int}(n_2/3)] \\
m_5 \in [0,n_1] \\
m_6 \in [0,n_2]
$$
where $\text{int}(n)$ stands for the integer part of $n$ (the "floor" function). By inverting my matrix, I obtain a solution generator:
$$
\begin{pmatrix}
m_1 \\
m_2 \\
m_3 \\
m_4 \\
m_5 \\
m_6
\end{pmatrix}
=
\begin{pmatrix}
1 & 1 & -2 & -3 & 1 & -1 \\
0 & -1 & 2 & 3 & -2 & 1 \\
0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
n_1 \\
n_2 - 2n_1 \\
m_3 \\
m_4 \\
m_5 \\
m_6
\end{pmatrix}
$$
Solutions are generated algorithmically for the fixed $n_i$ and varying $m_3$ to $m_6$ within their predetermined bounds (as given above; more efficient bounds could be determined).
All the valid solutions are those for which $m_i \in \mathcal{N}$, all other solutions (which in this case involve negative $m_i$) are discarded.