Find all natural (+0) solutions of a system of linear equations with natural (+0) coefficients 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.
 A: 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.
