Solving linear system of nonhomogenous equations that are known to have natural solutions
You can represent the system as a matrix equation as $Ax = b$, where $k^{th}$ row of $A$ is the coefficients $[a_{k1},a_{k2}...,a_{kn}]$, $x = [x_1, x_2,...,x_n]$ and $b = [b_1,b_2,...,b_m]$.
Unique Solution:
If the system has a unique solution, you should be able to get $A$ into row-echelon form. Then, you'll find a row with only one non-zero coefficient, corresponding to an atomic expression like $x_j = b_k$. Plug this in by dropping the $j^{th}$ column and the $k^{th}$ row and updating all other rows as $b_l' = b_l - a_{lj}b_k$. Recursively keep doing this until all values are found.
I finally remembered that the name of this procedure: Gaussian Elimination. Here is a nice concise description with a solved example.
Underdetermined System:
If there are multiple solutions, then the number of linearly independent rows is less than the number of variables. Since you seem to have many more equations than variables, this is great news. Using the row-echelon form, you can weed out redundant equalities and end up with a much smaller system of equations.
Now, let $$z_k = \min_{\substack{1 \leq i \leq m \\ a_{ki} \neq 0}}\left\{\frac{b_i}{a_{ki}}\right\},$$
so we get $1 \leq x_k \leq z_k$ as the feasible set for each variable. One way to shrink the brute-force approach is to order the variables by their upper bounds. If $z_{k_1} \geq z_{k_2} \geq ... \geq z_{k_n}$, suppose variables satisfy $x_{k_1} \geq x_{k_2} \geq ... \geq x_{k_n}$. Then, you can enumerate the solutions under this constraint by enumerating.
For example, let's say $z_1 = 3, z_2 = 2, ...$. Then, you would have $(3,2,..),...,(3,1,..),...,(2,2,..),...,(2,1,...),...,(1,1,...),...$.
With this ordering assumption you are searching over a much smaller set. Once you have that, you can check the permutations of these solutions and if they are feasible for all variables, they are valid solutions.
For example, for $x_1 + x_2 + x_3 = 3$, let's say we break the tie as $x_1 \geq x_2 \geq x_3$. Then, the constrained solutions are $(3,0,0), (2,1,0), (1,1,1)$ and the permutations of these give you the whole set since they are all the feasible due to $z_1=z_2=z_3$.