%Script file to produce the bifurcation diagram for
%  F(x) = x^2 + c using the orbit of 0, and values
%  of c:  -2 <= c <= 0.25

c=linspace(-2,0.25,1000);  %This is 1000 points for c
x=zeros(1,1000);           %An array of initial conditions
fprintf('Going into the initial orbit\n');

for i=1:100
   y=x.^2 + c;
   x=y;
end
fprintf('Going into new orbit\n');
%Now we're ready to store the results.
A=zeros(800,1000);
for i=1:800
   y = x.^2 + c;
   A(i,:)=y;
   x=y;
end

plot(A','k.','Markersize',4);