Are they using the same versions of SymPy?
e1 is quadratic only if f and e2 is linear in g and quadratic in f...it is a pretty simple system but there are lots of extraneous symbols. Even paring down the symbols to the system `(F1*(F2 + F3*f + f**2), G3*(G4 + G6*f**2 + G7*g + f*(G1 + G5*g - g))` takes a long time to solve and simplify.
That’s my point… BTW : with check=False, solving with Sympy’s solver is about as fast as with Maxima’s solver ; checking is slower, simplifying is extremely slow :
>>> from sympy import symbols, solve
>>> from time import time as stime
>>> f, q, r, h, i, g = map(symbols, 'fqrhig')
>>> e1 = -4*f**2*q - 5*f**2*r - 5*f*r + 4*q**2 + 3*q*r
>>> e2 = -4*f**2*h*q - 4*f**2*i*q - 8*f*g*q - 10*f*g*r - 8*f*h*q - 20*f*h*r - 10*f*i*r + 8*f*q**2/5 - 5*g*r + 4*h*q**2 - 5*h*r + 4*i*q**2 + 4*q**2/5 + 10*q*r + 5*r**2
>>> t0 = stime() ; Sol = solve([e1, e2], [f, g], check=False, dict=True) ; t1 = stime() ; print(t1 - t0)
0.6349234580993652
>>> t0 = stime() ; Chk = [[e.subs(s).simplify() for e in (e1, e2)] for s in Sol] ; t1 = stime() ; print (t1 - t0)
64.00839900970459
>>> Chk
[[0, 0], [0, 0]]
>>> t0 = stime() ; SSol = [{u:s[u].simplify() for u in s.keys()} for s in Sol] ; t1 = stime() ; print(t1 - t0)
1.9784221649169922
>>> SSol
[{f: -(5*r + sqrt(64*q**3 + 128*q**2*r + 60*q*r**2 + 25*r**2))/(8*q + 10*r), g: (-(4*q + 5*r)**2*sqrt(64*q**3 + 128*q**2*r + 60*q*r**2 + 25*r**2)*(160*h*q**3*r + 200*h*q**2*r**2 + 250*h*q*r**2 + 625*h*r**3 + 160*i*q**3*r + 200*i*q**2*r**2 + 250*i*q*r**2 + 625*i*r**3 + 64*q**4 + 880*q**3*r + 2400*q**2*r**2 + 2250*q*r**3 + 625*r**4) + (-80*h*q**2 - 250*h*q*r - 250*h*r**2 - 50*i*q*r - 125*i*r**2 + 16*q**3 + 20*q**2*r)*(1024*q**5 + 4608*q**4*r + 7680*q**3*r**2 + 5600*q**2*r**3 + 400*q**2*r**2 + 1500*q*r**4 + 1000*q*r**3 + 625*r**4))/(5*(4*q + 5*r)**2*(1024*q**5 + 4608*q**4*r + 7680*q**3*r**2 + 5600*q**2*r**3 + 400*q**2*r**2 + 1500*q*r**4 + 1000*q*r**3 + 625*r**4))}, {f: (-5*r + sqrt(64*q**3 + 128*q**2*r + 60*q*r**2 + 25*r**2))/(2*(4*q + 5*r)), g: ((4*q + 5*r)**2*sqrt(64*q**3 + 128*q**2*r + 60*q*r**2 + 25*r**2)*(160*h*q**3*r + 200*h*q**2*r**2 + 250*h*q*r**2 + 625*h*r**3 + 160*i*q**3*r + 200*i*q**2*r**2 + 250*i*q*r**2 + 625*i*r**3 + 64*q**4 + 880*q**3*r + 2400*q**2*r**2 + 2250*q*r**3 + 625*r**4) + (-80*h*q**2 - 250*h*q*r - 250*h*r**2 - 50*i*q*r - 125*i*r**2 + 16*q**3 + 20*q**2*r)*(1024*q**5 + 4608*q**4*r + 7680*q**3*r**2 + 5600*q**2*r**3 + 400*q**2*r**2 + 1500*q*r**4 + 1000*q*r**3 + 625*r**4))/(5*(4*q + 5*r)**2*(1024*q**5 + 4608*q**4*r + 7680*q**3*r**2 + 5600*q**2*r**3 + 400*q**2*r**2 + 1500*q*r**4 + 1000*q*r**3 + 625*r**4))}]
HTH,
“never returns” is a bit of an over statement. About half an hour is not eternity :
>>> from time import time as stime
>>> t0 = stime() ; Sol2 = solve([e1, e2], [f, g], dict=True) ; t1 = stime() ; print(t1 - t0)
1628.6711568832397