a = 3.56994567;
logistic = @(x,a) a.*x.*(1-x);

alpha = 2.8;
alpha = 2.502907875095892822283902873218215786381271376727149977336192056779235; 
%Feigenbaum alpha  from https://oeis.org/A006891

max_exponent = 12;

% Focus on the behaviour around 0.5...
% But you can increase the range up to [0 1]
x_values = linspace(0, 1, 1000000);
y_values = x_values;

last_n = 0;
for n=0:max_exponent
    n_iter = 2^n;

    for k=0:((n_iter-last_n)-1)
        y_values = logistic(y_values,a);
    end
    last_n = n_iter;

    plot(0.5+(x_values-0.5)*(-alpha)^(n), 0.5+(y_values-0.5)*(-alpha)^(n),'DisplayName',['-{\alpha}^{' num2str(n) '}f^{2^{' num2str(n), '}} (x/{\alpha}^{' num2str(n) '})']);
    hold on;
    xlabel('x')
    ylabel('Iterated maps')
    legend

    l_x = 5;
    xlim([0.5-l_x, 0.5+l_x]);
    ylim([-6, 4]);
    hold on;
    drawnow;
    pause(1)
end