In the following, a natural number is defined as a nonnegative integer.
Given a matrix $A$ consisting of only natural numbers, and a vector $\vec b$ of natural numbers, how do you find the set of solutions $\vec x$ of the equation $A \vec x = \vec b$ where all the components of $\vec x$ are natural numbers?
The approach described on Wikipedia for solving systems of Diophantine equations does not ensure that the solutions are nonnegative. Adding that constraint makes it a special case of integer linear programming, which is NP-Hard in general. Is there a polynomial time algorithm?
Indeed, this is a special case of the Knapsack Problem. It is NP-complete, so that there is no polynomial-time algorithm.
Reference: Theoretical Computer Science Stack Exchange.