% set up a collection of c values
x1 = linspace(-.75, 0.25, 7);
x2 = linspace(-1.25, -0.75, 22);
x3 = linspace(-1.4, -1.25, 22);
x4 = linspace(-1.75, -1.4, 22);
x5 = linspace(-1.75, -1.76, 22);
x6 = linspace(-2, -1.78, 22); 
% remove endpoints from the intervals of c values that are sampled
x1 = x1(2:end-1);
x2 = x2(2:end-1);
x3 = x3(2:end-1);
x4 = x4(2:end-1);
x5 = x5(2:end-1);
x6 = x6(2:end-1);

x = [x6 x5 x4 x3 x2 x1]; % gather collection of c values into one list
F = @(x,c) x^2+c; % our map of interest (the quadratic map)

n = 1000; % number of steps in each orbit
orbit = zeros(length(x),n+1); % matrix of orbits, each row corresponds to a c value and contains the orbit of 0 for that c-value

% compute the orbit for each c value stored in x
for i = 1:length(x)
    for j = 1:n
        orbit(i,j+1) = F(orbit(i,j), x(i));
    end
end

tailsize = 50; % we will plot the last 50 steps of each orbit
tailorbit = reshape(orbit(:, (n-tailsize+2):(n+1))', [1,length(x)*tailsize]); % some data wrangling to make the plot
scatter(repelem(x, tailsize), tailorbit, 10, 'filled'); % make a plot of points in the c-x plane
