Sympy reduction of a multivariate, bounded, under-determined system This is bad. I've been hacking away at Sympy for a few hours and the only solution I've come up with to yield the data I want is brute-force:
Code
from sympy import Interval, Naturals0, Q, Symbol, refine, solveset

N = 10
digit_interval = Interval.Ropen(0, N)
digits = digit_interval.intersect(Naturals0)  # equivalent to range(N)

# These don't seem to have any useful impact on solveset
assumptions = {
    'nonnegative': True,
    'noninteger': False,
    'real': True,
}

a = Symbol('a', **assumptions)
b = Symbol('b', **assumptions)
c = Symbol('c', **assumptions)
symbols = (a, b, c)

left = a / (b + c)
right = 3
soln = solveset(left - right, symbols[-1], digits)

# Useless - has no effect
for sym in symbols:
    soln = refine(soln, Q.positive(N - sym))

print(f'Set solution to {left}={right} for {symbols[-1]} is {soln}')


# Number of variables is so small that recursion is fine(ish)
def recurse_determine(inner_soln, vars=(), depth=0):
    if depth == len(symbols) - 1:
        yield from (vars + (v,) for v in inner_soln)
    else:
        for val in digits:
            yield from recurse_determine(
                inner_soln.subs(symbols[depth], val),
                vars + (val,),
                depth+1,
            )


for values in recurse_determine(soln):
    print(', '.join(str(v) for v in values))

The crux of the problem is that it seems difficult (impossible?) to have Sympy find a solution to such a multivariate integer problem and meaningfully constrain it to known bounds, so I'm left doing the brute-force approach: iterate through every possible value of every variable, and only when a (probably invalid) selection of variable values has been substituted can I ask Sympy to tell me that the result is an empty set. It works (for about the worst definition of &quot;works&quot;) but there has to be a better way.
I should add that the bounds will remain the same but the equation will be subject to change in the future, so a suggested solution that replaces Sympy with (e.g.) hard-coded Numpy or bare Python will not help me.
Output
Set solution to a/(b + c)=3 for c is Intersection(FiniteSet(a/3 - b), Range(0, 10, 1))
0, 0, 0
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
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.