Given an instance of 3-sat, create a variable (in the column vector x) for each variable and its negation in the boolean expression. To illustrate, let z and ~z be two such variables. Perhaps these are represented by x7 and x8 in the column vector x. Create four inequalities, four rows in the matrix A, that force z ≥ 0, ~z ≥ 0, z + ~z ≥ 1, and - z-~z ≥ -1. Now z can take on any value in [0,1], with ~z = 1-z. Yet all entries are integers, hence z is 0 or 1, and ~z is 1-z.
If a clause looks like x|y|~z, add another row to the matrix that forces x+y+~z ≥ 1. This system of inequalities has an integer solution iff the boolean expression is satisfied.
This problem is atypical, since it is easily shown to be np hard, while the task of proving it is in np (which is usually so straightforward that we don't even talk about it) is more difficult. Given a matrix A and a column vector b, an ntm could try millions of integer values for x1 through xn concurrently, but it still might not find a solution in polynomial time. It could try larger and larger values of x, but where does it end? Is the problem even recursive?
If we could put a bound on x, and if the bound isn't ridiculously large, the ntm can try all values of x concurrently, and its lucky guess will be rewarded, making this an np problem. How can we constrain the integers that make up the column vector x?
Each inequality, represented by a row in the matrix A, defines a hyperplane in n space. The values of x are forced to lie on or above this hyperplane. Together, the m hyperplanes define a polyhedron that encloses a volume, that contains x. We need to establish the boundary of this polyhedron in n space.
At this point I defer to an area of mathematics called linear programming, or linear optimization, which explores the interior of these shapes, and searches for extremal points within these confines. Polynomial time procedures determine whether the region is closed, and set limits on the coordinates. These procedures involve linear transformations and matrix operations. Integers can become rational numbers, as hyperplanes intersect, but cramer's rules put a bound on the common denominator, according to the determinant of the matrix. When written in binary, the length of the number that is the determinant can be the square of the length of the largest number in the matrix. This leads to the bound on x, and it's no worse than the square of the problem instance. Therefore all values of x can be tried in parallel, and the problem is np, and np complete.
As you can see, I've skipped past a lot of details. Someday I'll bring linear optimization online, and when I do, I'll refer you to the appropriate theorems. In the meantime, this outline will have to suffice.
Slip in an equation if you like; it's the same as two inequalities: Ai.x ≤ bi and Ai.x ≥ bi.
Strict inequality doesn't change anything either. If Ai.x > bi, rewrite it as Ai.x ≥ bi+1. Thus all flavors of simultaneous equations and inequalities are permitted, and the task of finding integer solutions remains np complete. Of course, if only equations are involved, a rational solution space is quickly found, and it's easy to check whether this includes an integer solution.