% Choose the number of c's and number of iterations and start value of x.
Nx=100;
Nc=20000;
c=linspace(0.25,-1.5,Nc)';
x=zeros(size(c,1),Nx);
x(:,1)=0*ones(size(c,2),1);

%Iterate the function x^2+c
for k=2:Nx
    x(:,k)=x(:,k-1).^2+c;
end

% Plot the last 20 iterations
Nitplot=20;
plot(c,x(:,end-Nitplot:end),'+','markersize',2);


% Label axes and add a legend:
xlabel('c'); ylabel('F^N'); title('Orbits of F=x^2+c');
h=get(gca,'children');
legend(h(1:11),num2str(Nx),num2str(Nx-1),num2str(Nx-2),num2str(Nx-3),...
    num2str(Nx-4),num2str(Nx-5),num2str(Nx-6),num2str(Nx-7),...
    num2str(Nx-8),num2str(Nx-9),num2str(Nx-10),1);