In sympy if you want to find integer solutions to equations then you should use diophantine. It doesn't handle systems of equations but you can put the solution from one equation into the other and call diophantine again:
In [69]: eq1 = x + y + z - 100
In [70]: eq2 = 10*x + 5*y/2 + z/2 - 100
In [71]: sol = diophantine(eq1, t, syms=[x, y, z])
In [72]: sol
Out[72]: {(t₀, t₀ + t₁, -2⋅t₀ - t₁ + 100)}
In [73]: [xt, yt, zt], = sol
In [74]: eq3 = eq2.subs({x:xt, y:yt, z:zt})
In [75]: eq3
Out[75]:
23⋅t₀
───── + 2⋅t₁ - 50
2
In [76]: t1, t2 = eq3.free_symbols
In [77]: [t1s, t2s], = diophantine(eq3, z, syms=[t1, t2])
In [78]: rep = {t1:t1s, t2:t2s}
In [79]: (xt.subs(rep), yt.subs(rep), zt.subs(rep))
Out[79]: (4⋅z₀ - 100, 500 - 19⋅z₀, 15⋅z₀ - 300)
The solution here is in terms of an integer parameter z0. This gives the set of solutions to the two equations but you also have the requirement that x, y, z are positive which constrains the possible values of z0:
In [80]: ineqs = [s.subs(rep) > 0 for s in [xt, yt, zt]]
In [81]: ineqs
Out[81]: [4⋅z₀ - 100 > 0, 500 - 19⋅z₀ > 0, 15⋅z₀ - 300 > 0]
In [82]: solve(ineqs)
Out[82]:
500
25 < z₀ ∧ z₀ < ───
19
In [83]: 500/19
Out[83]: 26.31578947368421
We see that z needs to be 26 which gives a unique solution for x, y and z:
In [84]: z, = ineqs[0].free_symbols
In [85]: (xt.subs(rep).subs(z, 26), yt.subs(rep).subs(z, 26), zt.subs(rep).subs(z, 26))
Out[85]: (4, 6, 90)