Sympy reduction of a multivariate, bounded, under-determined system
You're more running into the limitations of SymPy than anything. Most of the things you are doing should work, they just aren't implemented presently.
Your best bet for solving this kind of problem in a way that is smarter than naive brute force iteration of all possible values is to use the Diophantine solvers:
&gt;&gt;&gt; diophantine(a/(b + c) - 3, t, syms=(a, b, c))
{(3*t_0, t_1, t_0 - t_1)}

This tells you that a solution set is
a = 3*t0
b = t1
c = t0 - t1

where t0 and t1 are arbitrary integers.  Given the form of the solution, it is quite easy to see that for a, b, c to be in the range [0, N], t1 must be in [0, N] and t0 must be in [0, N/3].
&gt;&gt;&gt; N = 10
&gt;&gt;&gt; for t0 in range(0, N//3 + 1):
...     for t1 in range(0, N + 1):
...         a, b, c = 3*t0, t1, t0 - t1
...         if all(0 &lt;= i &lt;= N for i in (a, b, c)) and b + c != 0:
...             print(a, b, c)
3 0 1
3 1 0
6 0 2
6 1 1
6 2 0
9 0 3
9 1 2
9 2 1
9 3 0

Actually in this case, brute force is also fine because there are only N^3 = 1000 combinations, but Diophantine reduced it to only N^2, so that's an improvement asymptotically. There's no need to do it recursively. You can build a substitution dictionary and pass it to subs (like expr.subs({a: 1, b: 2, c: 3})).
You can see the documentation for more information on what kinds of Diophantine equations are implemented. I'll leave it to you to implement the general code to do the above. It really depends on what form your equations take in general as to what their Diophantine solutions will look like, including whether SymPy can solve them at all.