I have three huge degree 31 bivariate polynomials (20,000 characters long each) I want to plot, but I keep getting a lot of noise in my plot. I can't upload it, but the point is that in some regions I just get colorful noise. I've tried defining the polynomials over RealField(n) and increasing the number of plot_points, but neither of these approaches work. Any ideas on how to work around this? Thanks.

**Edit:** Tried using sympy's plot_implicit and it's so (SO!) slow. Then used numpy's contour_plot and it's fast, but has the same problem as sage.

Here's the code that produces the polynomials and plot. Be patient as it could be a bit slow (depending on your machine).

```
plot = Graphics()
m = 32-1
D = [(i,j) for i in range(0,m+1) for j in range(0,m+1) if i+j<m+1]
#Polygon
P = Polyhedron( vertices= [(0, 0), (32, 44), (23, 0), (10, 14), (2, 3)] )
points = P.integral_points()
plot_pts = point(points, rgbcolor=(0, 0, 0), size = 20).plot()
plot_np = P.plot(fill = False, point=False, line='black')
M = matrix(ZZ, len(points), len(D), 0)
for row_num, row in enumerate(points):
for col_num, column in enumerate(D):
i, j = row
a, b= column
#Matrix for interpolation:
M[row_num, col_num] = (i^a)*(j^b)
R = PolynomialRing(QQ, 2, 'xy')
S = PolynomialRing(RealField(500), 2, 'uv')
x, y = R.gens()
u, v = S.gens()
K = M.right_kernel()
Kdim = K.dimension()
print(Kdim)
if Kdim > 0:
for l in range(Kdim):
K_basis = K.basis()[l]
#Writing the interpolating polynomial
f=0
for order, bidegree in enumerate(D):
a, b = bidegree
f += list(K_basis)[order] * u^a * v^b
F = f.factor()
f = F[6][0]
cols = ['red', 'blue', 'green']
interpolation = implicit_plot(f, (u,-1,34), (v,-1,45), plot_points=100, color=cols[l])
plot += interpolation
plot += plot_pts + plot_np
plot.show(figsize=10)
```

**Edit 2:** Using the mpmath library in Python and with the aid of Sébastien's code below I wrote a routine that allows us to control the root finding method and precision of our computations. I tried several methods, secant (default), newton, hailley, mnewton, etc. without success. Changing the precision and tolerance of the root finding function from mpmath didn't help either. I think this polynomial just behaves too wildly in the region of the plot.

Here's the code and relevant documentation for the "mpmath.findroot" function:

http://mpmath.org/doc/current/calculus/optimization.html :

```
from sympy import *
import matplotlib.pyplot as plt
import mpmath as mp
mp.dps = 100
stepx = 0.1
stepy = 0.5
xrange = np.arange(7.5,12.5, stepx)
yrange = np.arange(0, 5, stepy)
def plot_roots_of_f(f):
L = []
for u in xrange:
for v in yrange:
Root = mp.findroot(f.eval({x:u}),v, method = 'mnewton', tol = E-60, verbose=False,verify=False)
L.append([u,Root])
return plt.plot(*zip(*L),linestyle='None', marker='.')
plot_roots_of_f(f)
plt.show()
```