function M = mandelbrot( d, it )

 re = [-2.0 1.0]; 				# Domain
 im = [-1.5 1.5];  				# Range

## d = 0.05;					# Step size
## it = 16;               			# Iteration depth

  m_re = (re(2)-re(1))/d;
  m_im = (im(2)-im(1))/d;

 ## Prepare the complex plane
 regrid(m_re)=0;
 imgrid(m_im)=0;

 for j=1:m_re
   regrid(j) = re(1) + (j-1)*(re(2)-re(1))/(m_re-1) ;
 end
 for j=1:m_im
   imgrid(j) = im(1) + (j-1)*(im(2)-im(1))/(m_im-1) ;
 end

 ## Generate Mandelbrot set M
 M(m_re,m_im)=0;
 c_re=0;
 c_im=0;

 for j=1:m_re
  for k=1:m_im
    c_re = regrid(j);
    c_im = imgrid(k);

    re_xn = 0 ;
    im_xn = 0 ;
    re_xn1 = 0 ;
    im_xn1 = 0 ;
    n = 0 ;
    
    while sqrt(re_xn*re_xn+im_xn*im_xn)<=2 && n<it
      re_xn1 = re_xn*re_xn - im_xn*im_xn + c_re ;
      im_xn1 = 2*re_xn*im_xn + c_im ;
      re_xn = re_xn1 ;
      im_xn = im_xn1 ;
      n++ ;
    end

    M(j,k) = n ;

  end
 end