Let a chemical system be defined as $${A<=>B<=>C}$$
Then the stoichiometry is given as $$S=\begin{bmatrix} -1& 1& 0& 0\\ 1& -1& -1 & 1\\ 0 & 0 & 1 &-1\end{bmatrix}$$
Then the conservation law can be found as (MATLAB) $$C=rref(null(S′,′r′)′)′$$ Then $$C=\begin{bmatrix} 1\\1\\1\end{bmatrix}$$ For a closed chemical system, if the initial state, $\vec{X}(0)$ ($\vec X(t)$= number of molecules A, B, C at time t), is given, then for all time, states will satisfy the conservation. That means $$C^T\vec{X}(t)=C^T\vec{X}(0)\hspace{1cm} \forall t$$
In the case of chemical system explained here, for an initial state $\vec X(0)=\begin{bmatrix} 2&1&0\end{bmatrix}^T$, $$\begin{bmatrix} 1&1&1\end{bmatrix}\vec{X}(t)=\begin{bmatrix} 1&1&1\end{bmatrix}\begin{bmatrix} 2\\1\\0\end{bmatrix}=3$$
Then the question is how can I identify all possible $\vec X(t)$ such that $x_i\in \mathbb Z_+ \ge 0$
What I need is a general strategy that can find all columns of $\{X\}$ given below.
$$\{X\}=\begin{bmatrix} 0 & 1& 2& 3& 0 & 0 & 0 & 1 & 1 & 2\\ 0 & 0& 0& 0& 1 & 2 & 3 & 1 & 2 & 1\\ 3 & 2 & 1 & 0 & 2 & 1 & 0 & 1 & 0 & 0 \end{bmatrix}$$
As you can see $$C^T\{X\}=\begin{bmatrix}3&3&3&3&3&3&3&3&3&3\end{bmatrix}$$
There is no other positive integer combination is possible other than these columns. So basically I need to solve an equation
$$Ax=b\hspace{1cm} A\in \mathbb Z_+^{d\times n}\hspace{1cm} x\in \mathbb Z_+^{1\times n}\ge 0 \hspace{1cm} b\in \mathbb Z_+^{1\times d} \ge 0\hspace{1cm} d< n$$
A similar question: Getting all integer solution to undertermined linear system
Finding all the solutions to a given system of $\mathbb{Z}$-linear equations is a hard mathematical problem. In your application, is $A$ always $\begin{bmatrix}1&1&1\end{bmatrix}$? If so, then for an initial state $X\in\mathbb{Z}_+^3$, the solutions to your system are precisely the partitions of $X_1+X_2+X_3\in\mathbb{N}$ into $3$ parts.
For general matrixes $A$ and right-hand sides $b$, you can use the software LattE which counts the number of solutions.