mandel(x1,x2,y1,y2,r,n,s,it)

function mandel(x1,x2,y1,y2,r,n,s,it)

% MANDEL    Draws the Mandelbrot set.
%   MANDEL(x1,x2,y1,y2,r,n,s,it) draws the Mandel-
%   brot set in the region R = [x1,x2] x i [y1,y2]. 
%   The number n defines a grid in R with n^2 points. 
%   The number it denotes the number of iterations 
%               z -> z^r + c
%   taken by MANDEL.
%   The result is better for large values of n and it.
%
%   MANDEL does not draw each iteration, but only each
%   s-th iteration. A good choice is 5 <= s <= 10.
%   
%   You may start with x1 = -2, x2 = 1, y1 = -1.2, 
%   y2 = 1.2 and then "zoom in" by just specifing
%   the new region with the cursor in the current 
%   figure.
%   
% ==========================================================================
% intialization
tic,mandelcalc(x1,x2,y1,y2,r,n,s,it),toc

% refine 
refine = input('\n Do want to zoom in ? Type 1 for yes, 0 for no ... ');
while refine == 1
    [X,Y] = ginput(2);
    n = input('\n Number of grid points ... ');
    it = input('\n Number of iterations ... '); 
    mandelcalc(X(1),X(2),Y(1),Y(2),r,n,s,it)
    refine = input('\n Do want to zoom in ? Type 1 for yes, 0 for no ... ');
end

% ==========================================================================
function mandelcalc(x1,x2,y1,y2,r,n,s,it)

[X,Y] = meshgrid(0:n-1,0:n-1);
c = x1+(x2-x1)/n*X+i*(y1+(y2-y1)/n*Y);
c = reshape(c,n^2,1);
Z = 0*c;

indices = find(abs(Z)<2);
p = 0;
figure
while (~isempty(indices)) & (p < it)
    indices1 = indices;
    for k = 1:s
        Z(indices) = Z(indices).^r+c(indices);
    end
    indices = find(abs(Z)<2);
    plot(real(c(setdiff(indices1,indices))),imag(c(setdiff(indices1,indices))),...
        '.','Color',[((it-p)/it)^2 ((it-p)/it) (it-p)/it],'MarkerSize',0.5)
    hold on
    axis('off'), axis([x1 x2 y1 y2]), axis('equal'), pause(0.01)
    p = p+s;
end