aを定数とすると、与式が自然数ならば
b+a^2+a≧ab^2+b+7
⇔ab^2-a^2-a+7≦0
⇔-√(a+1-7/a)≦b≦√(a+1-7/a) (∵a>0)
bは自然数だから
1≦b≦√(a+1-7/a)
特にa≦10000のときb≦100<100.0…
# Python 3
for a in range(1,10001):
for b in range(1,101):
k=(b+a*a+a)/(a*b*b+b+7)
if k==int(k):
print(a,b,k)
print('done')
数秒で次の出力を得た
11 1 7.0
17 2 4.0
27 3 3.0
49 1 43.0
done
これ以外に解はないと考えられる