c = 1.5i;
Finv_c = @(x) sqrt(x-c);

x0 = 0;
n = 10000;

orbit = zeros(1,n+1);
orbit(1) = x0;
for i = 1:n
    orbit(i+1) = (2*randi([0 1])-1)*Finv_c(orbit(i));
end

plot(orbit(101:n), '.')