c = -1; % 
F_c = @(x) x^2 + c;

n = 1000; % grid size
N = 100; % max iterations

x = linspace(-2,2,n);
[a,b] = meshgrid(x);
z = complex(a,b);
gridcolor = zeros(n,n);

for j = 1:n
    for k = 1:n
        orbit = zeros(1, N+1);
        orbit(1) = z(j,k);
        for i = 1:N
            if abs(orbit(i)) > max(abs(c), 2)
                gridcolor(j,k) = 1;
                break
            end
            orbit(i+1) = F_c(orbit(i));
        end
    end
end


scatter(reshape(real(z), [n*n, 1]), reshape(imag(z), [n*n, 1]), 10, reshape(gridcolor, [n*n, 1]), 'Filled');

