p=20 # experimente usar outros valores para p e b
b=20
Z=0*(1:(p+1))
Z[1]=0
c=.10+.3*1i  # experimente usar outros valores aqui
h=b/p; Y=0

f<-function(U){p=U^2-U+c;p } 


for (i in 1:p){ Y=Y+h*f(Y); Z[i+1]=Y}

print("Valores de Z_(n+1)=Z_n^2+c, iniciando como Z_0=0, para c=.10+.3*1i."); Z

plot(Re(Z),Im(Z),'l',ylim=c(0,0.4),col="blue")
legend("topright",legend="Tragetória de Z_(n+1)=Z_n^2+c, iniciando como Z_0=0, para c=.10+.3*1i.", lty=1, col="blue", lwd=1, bty="n")

# Mandelbrot

s=seq(-2,1,by=0.01);  t=seq(-1.5,1.5,by=0.01); n=length(s)
A=matrix(0,n,n)

for (i in 1:n){ for (j in 1:n){ Y=0; c=s[i]+t[j]*1i

f<-function(U){p=U^2-U+c;p } 

for (l in 1:p){ K=abs(Y); Y=Y+h*f(Y)}
A[i,j]=K
}}

print("Região que contém os valores de em que Z_(n+1)=Z_n^2+c é limitada."); image(A,axes=FALSE)