What algorithms exist for solving natural number linear systems? I'm looking at the following problem:

Given $n$-dimensional vectors of natural numbers $v_1, \ldots, v_m$ and some input vector $u$, is $u$ a linear combination of the $v_i$'s with natural number coefficients?

i.e. are there some $t_1, \ldots, t_m \in \mathbb{N}$ where $u = t_1 v_1 + \dots + t_m v_m$?

Obviously the real-number version of this problem can be solved using Gaussian elimination. I'm wondering, has the integer version of this problem been studied? What algorithms exist to solve it?

Note that this is using natural numbers, but not modular arithmetic, so this is somewhat separate from the Chinese Remainder Theorem and systems like that. Also, it seems related to Diophantine equations, but I'm wondering what has been done in the case where only non-negative integers are considered?  This is also reminiscent of a multi-dimensional subset-sum problem, generalized to allow us to take an arbitrary number of copies of each vector.  It also seems related to testing whether $u$ is an element of the lattice generated by $v_1,\dots,v_m$, except that here we only allow linear combinations with non-negative coefficients.

For anyone interested, this is motivated by looking at whether a Parikh vector is in a linear set, as in Parikh's Theorem.

In particular, I'm interested in an algorithm that could solve the problem using only natural number operations, avoiding going into the reals/floating point numbers.
Your problem is NP-complete, by reduction from Subset Sum (it is in NP since the fact that everything is non-negative bounds the coefficients of the solution sufficiently well). Given an instance $S = \{s_1,\ldots,s_n\}, T$ of Subset Sum (is there a subset of $S$ summing to $T$?), we construct an instance $v_1,\ldots,v_{2n},u$ of your problem as follows. For each $1 \leq i \leq n$, we put $v_i$ to be the vector with two non-zero entries: $v_{i,i} = 1$ and $v_{i,n+1} = s_i$, and $v_{n+i}$ to be the vector with a unique non-zero entry $v_{n+i,i} = 1$. The target vector is $u=1,\ldots,1,T$. Each natural combination of $v_1,\ldots,v_{2n}$ equal to $1,\ldots,1,\ast$ must select exactly one of each of $v_i,v_{n+i}$, and so encodes a subset of $S$ whose sum is the value of the last component.