function W=mandelprot1(X,Y,a,n_iter)
C=X+i*Y; W=zeros(size(X)); Z=1;
for k=1:n_iter, 
    Z=(2*Z.^3+C)./(3*Z.^2); 
    i0=find(abs(Z)>a); W(i0)=k; Z(i0)=NaN;
end
i0=find(W==0); W(i0)=NaN;