Consider the following equation:
$$ \begin{pmatrix} n_1 \\ n_2 \end{pmatrix} = \begin{pmatrix} 1 & 1 & 0 & 0 & 1 & 0 \\ 2 & 1 & 2 & 3 & 0 & 1 \end{pmatrix} \begin{pmatrix} m_1 \\ m_2 \\ m_3 \\ m_4 \\ m_5 \\ m_6 \end{pmatrix} $$
where $n_i \in \mathcal{N}$ are known and $m_i$ are unknown. I would like to find all the possible solutions $m_i \in \mathcal{N}$. Note the trivial solutions $m_i = n_i = 0$ and $m_{1,2,3,4} = 0; m_5 = n_1; m_6 = n_2$. Also note I have considered $\mathcal{N}$ to contain zero.
What's an algorithmic way to get all possible solutions compatible with given $n_i$? Something implementable in, e.g., Python, is OK for my purposes.
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.