% Plot the Bifurcation diagram of the logistic map

mMin    = 0;    % Smallest value of m for which to plot
mMax    = 2;    % Largest value (+step size) of m for which to plot 
steps   = 200;  % Number of subdivisions on the m-axis
preiter = 100;  % Number of pre-iterations (not plotted)
iter    = 30;   % Number of iterations to plot

h = (mMax - mMin)/steps; % Step size on the m-axis

% Pre-allocate arrays for speed only
xx = zeros(1,steps*iter); 
yy = zeros(1,steps*iter); 
y  = zeros(1,iter);

for j=0:(steps-1)
  m = j*h;
  x = m*ones(1,iter);
  y = x; % prealloate for speed;

  % First, pre-iterate the logistic map without storing values 
  y(1) = 1;
  for k=1:preiter
    y(1) = 1 - m*y(1)^2;
  end

  % Now, store iter more iterations
  for k=2:iter
    y(k) = 1 - m*y(k-1)^2;
  end

  % Put the coordinates in the coordinate vectors for plotting
  xx((j*steps+1):(j*steps+iter)) = m*ones(1,iter);
  yy((j*steps+1):(j*steps+iter)) = y;
end

plot(xx,yy,'.')