function mandelpic

[X,Y] = meshgrid(-3:.001:1,-2:.001:2);
Z=X+1i*Y;

[p,dp]=mandelbrot(Z,7);

contour(X,Y,abs(Z.*dp./p),2.^[-4:1:8],'k');
axis('square');
axis image, axis off;

end

function [p,dp]=mandelbrot(z,n)
  if n<=0
    p=ones(size(z));
    dp=zeros(size(z));
  else
    [p0,dp0]=mandelbrot(z,n-1);
    p=z.*p0.^2+1;
    dp=p0.^2+2.*z.*p0.*dp0;
  end
end