using Plots
n = 5000
nit = 60
L = 3
x = LinRange(-L/2,L/2,n)
y = LinRange(-L/2,L/2,n)
X = repeat(x', n,1)
Y = repeat(y, 1,n)
Z = X .+ Y.*im
for it in 1:nit
@. Z = Z - (Z^7 - 1)/(7*Z^6)
end
phi = [angle(Z[i,j]) for i in 1:n, j in 1:n]
heatmap(x,y,phi, fill=true, c=:magma, size=(1000,1000), legend=:none, ticks=false)