import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import quad
from scipy.interpolate import interp1d
from scipy.optimize import curve_fit, fsolve
n = 4000 # number of points across plane
nit = 20 # number of Newton iterations
x = np.linspace(-2,2,n) # x points
y = np.linspace(-2,2,n) # y points
X,Y = np.meshgrid(x,y) # make 2-D array of x,y pts for plotting
Z = X + 1j*Y # 2D array of initial guesses
for it in range(nit) : # do Newton iterations (on the whole array at once)
Z = Z - (Z**4 - 1)/(4*Z**3)
plt.figure(figsize=(10,10))
phi = np.angle(Z)
plt.contourf(X,Y,phi,cmap='jet');