Vectorizing a nested loop algorithm which produces the Mandelbrot Set

I have a working piece of code that produces a rough plot of the Mandelbrot set (shown below) that I am trying to vectorize (also shown below). My indexing is flawed I believe in the vectorized code since the picture looks nothing like the Mandelbrot set. I would appreciate it if someone could point me to the problem in my vectorized code.

The basic version creates this:

The vectorized version creates this:

Please help.

%%%%%%%%%%%%THE WORKING NESTED LOOP CODE - BEGIN %%%%%%%%%%%%%%
% We want to iterate f(z) = z^2 + r

figure; hold on;
NumIts = 20
a = -2;
b = -2; 
c = 2; 
d = 2;
M = 100; %square root of the number of points
R = 10;

C = colormap(hot(NumIts));
% Nested loop generates a lattice (p,q)
ii=1;
for p = 1:M
  for q = 1:M
    % pick a point in C
    k = a+(c-a).*p./M;
    l = b+(d-b).*q./M;

    % specify the intial point, 0 on the orbit 
    x = 0;
    y = 0;

    % Compute at most NumIts iterations on the orbit of 0
    n = 1;
    while n < NumIts
      %FB{ii} = [x(Ix) y(Ix)]l;
      newx =  x.*x - y.*y - k; %This is the real part of z^2 - lambda
      newy =  2.*x.*y - l;    %The imaginary part

      x = newx;
      y = newy;
    
      
      mm(n) = x.*x + y.*y;
      myinfo{n} = [x y l k mm(n) n];
      B{ii}(n,:) = [x y];

      % If the iterate has escaped color it
      if x.*x + y.*y > R
        plot(p,q,'o','MarkerSize',5,'Color',C(n,:));
        n = NumIts;
      end
      n = n+1;
      ii =ii+1;
    end
  end
  disp(M-p)
end

%%%%%%%%%%%%%%WORKING NESTED LOOP CODE END %%%%%%%%%%%%%%%%%%%%%%%

Below is the vectorized code

%%%%%%%%%%%%%%Vectorized code that should produce a picture of the Mandelbrot set %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   % We want to iterate f(z) = z^2 + r

  % Set Constants
  clear

  NumIts = 20; % Number of iterations
  % Constants for the labmda value
  a = -2;
  b = -2; 
  c = 2; 
  d = 2;
  M = 100; % square root of the number of points to iterate
  R = 10; % the diameter of divergence

  % C is 3xNumIts matrix of RGB values
  C = colormap(hot(NumIts));
  figure; hold on;

  % To vectorize a calculation on a grid of points
  % Create a 3-d lattice of (p,q) pairs
  p = [1:M]';
  P = repmat(p,1,M);

  q = 1:M;
  for i = 0:M-1
    v = q+i;
    f = find(v>M);
    v(f) = mod(v(f),M);
    Q(:,i+1) = v;
  end

  % pick a point in C (he says P)
  k = a+(c-a).*P./M;
  l = b+(d-b).*Q./M;

  % specify the intial points 0 on the orbit 
  x = zeros(M);
  y = zeros(M);
  
  % Compute at most NumIts iterations on the orbit of 0
  n = 1;
  ItDist    = zeros(M);
  ItDistColor = zeros(M,M,3);
  Ix = [1:M^2]';
  Out = zeros(size(Ix));
  PointNorm = zeros(size(Ix));  
    
  % An array of input (x,y,l,k), output (PointNorm),
  % index and classification columns
  A = [x(Ix) y(Ix) l(Ix) k(Ix) PointNorm(Ix) Ix Out];

  while n < NumIts
    
    % Iterate this function z -> z^ + (k,l)
    newx = A(:,1).*A(:,1)-A(:,2).*A(:,2)-A(:,4);
    newy = 2.*A(:,1).*A(:,2) - A(:,3);

    % Update the iterate coordinates
    A(:,1) = newx;
    A(:,2) = newy;

    % Calculate the norm of the image of the iteration
    A(:,5) = A(:,1).*A(:,1)+A(:,2).*A(:,2);
    
    % Update the point classification variable to 1
    K = find(A(:,5) > R); 
    A(K,7) = 1;

    % Grab the large-distance coords 
    [I,J] = ind2sub([M,M],A(K,6));

    % Check that they're correct
    if find(A(find(A(:,7)),5)<R)
      disp('I & J indices not right');
      return;
    end

    % Plot them in a color matching the iteration depth
    ColorVec(1,1,:) = C(n,:);
    ItDistColor(I,J,:) = repmat(ColorVec,length(I),length(J));

    % Reduce to the lattice points whose image remains inside the
    % diameter of divergence
    A = A(find(A(:,5)<=R),:);
    
    % Debugging device
    BB{n} = A;

    n = n + 1        
  end

  % Plotting
  for i = 1:M; for j = 1:M;
    plot(i,j,'o','MarkerSize',5,'MarkerFaceColor',ItDistColor(i,j,1:3));
  end; end;