The question seems to be about programs that can compute the nullspace of an integer matrix exactly. (Though there seems to be the constraint that all the entries in the solution vector are positive, which I'm going to ignore.) The best software you might want to use could depend on the size of your matrices, which is a function of: (1) the parameters $M$, $N$, (2) the number of bits of the entries, and (3) the sparseness of the matrix.
IML (Integer Matrix Library) -- this is a free C library, but what it does is compute the rational kernel of a system defined over the integers, so it isn't directly useful for your problem. Anyway, in my benchmarks, this is the fastest library available for computing the rational kernel if your matrix is large (both $M$ and $N$), dense, and the entries are also large. When entry sizes get to several hundred bits or larger, IML seemed to me in 2007 to be by far the fastest library in existence (easily beating Mathematica, Magma, NTL, etc.). IML is used part of Sage and used for the kernel function. Getting going with IML standalone can be hard (though I contributed code a few years ago to make it easier), but if you just want a library to link into an existing C program, you'll definitely want to look into IML.
PARI -- this free calculator (which is also usable as a C library) is extremely good for computing kernels of small matrices, e.g., when the entries aren't too big and the number of rows and columns aren't too large either (say, less than 100). It's optimized for this small case, since that's what comes up in class group computations, among other things. The command to use is matkerint, e.g., "matkerint([1,2,3;4,5,6;7,8,9])".
SAGE -- this free software uses some combination of IML, PARI, Linbox, NTL and a bunch of new code to compute kernels of integers matrices. In some cases (e.g., nearly square matrices) it is very good, and it could be good in far more cases but I got lazy and didn't implement fast algorithms for all shapes of matrices. Type "A = matrix(ZZ,3,3,[1..9])" to make an example integer matrix; then "A.kernel()" will find the kernel over ZZ (by default, contrary to what another answer on this page claims). If you type "A._k[press tab key]" you'll see that you can directly call two other routines for computing kernels, one involving PARI, and the other a $p$-adic algorithm. Behind the scenes, really the Hermite normal form of the matrix is being computed using the algorithm described here. The HNF algorithm in Sage makes fundamental use of IML in several ways behind the scenes.
Magma -- non-free and closed source, but has a fast tuned Hermite normal form algorithm presumably similar to the one in Sage, which for certain (not all) ranges of inputs is the best available.
Mathematica -- in my benchmarks (from 2007), as the matrices get very big this gets totally useless compared to Magma or Sage, but it is fine for small matrices. Similar remarks apply about Maple, at least as of 2007.
Here's a potentially helpful example in Sage, which illustrates some of the computations discussed here and elsewhere, and explains how to work with some matrices and modules in Sage:
sage: A = random_matrix(ZZ,5,2); B = random_matrix(ZZ,5,2)
sage: A.kernel()
Free module of degree 5 and rank 3 over Integer Ring
Echelon basis matrix:
[ 1 60 0 1 30]
[ 0 66 0 1 33]
[ 0 0 1 0 0]
sage: A.kernel().intersection(B.kernel())
Free module of degree 5 and rank 1 over Integer Ring
Echelon basis matrix:
[ 19 -378 -93 -4 -189]
sage: VA = A.change_ring(QQ).kernel(); VB = B.change_ring(QQ).kernel()
sage: W = VA.intersection(VB); W
Vector space of degree 5 and dimension 1 over Rational Field
Basis matrix:
[ 1 -378/19 -93/19 -4/19 -189/19]
sage: W.intersection(ZZ^5)
Free module of degree 5 and rank 1 over Integer Ring
Echelon basis matrix:
[ 19 -378 -93 -4 -189]