Producing a bifurcation diagram -- How to make this code faster?

So the code produces a bifurcation diagram, however its taking really long to run. I was thinking if the same x value is produced on each iteration then the iteration needs to stop instead of being plotted on top of the old one. Unsure how you would code this though?

%Final plot
%Parameters -- Keep as constants for now.
d = 0.005;
%Parameters - Keep values of b and r constant.
b=8/3;
r=28;
for a = 0.2:0.001:1
%Initial conditions
xold=1.2;
yold=1.3;
zold=1.6;

for n=1:10000000
    xnew=xold+((d*a)*(xold-yold));
    ynew=yold+(d*((r*xold)-(zold*xold)-yold));
    znew=zold+(d*((xold*yold)-(b*zold)));
    xold=xnew;
    yold=ynew;
    zold=znew;
end

%Set as though xnew, ynew and znew are initial conditions
x(1)=xnew;
y(1)=ynew;
z(1)=znew;

for n=1:100
    x(n+1)=x(n)+(d*a)*(x(n)-y(n));
    y(n+1)=y(n)+(d)*((r*x(n))-(z(n)*x(n))-y(n));
    z(n+1)=z(n)+(d*((x(n)*y(n))-(b*z(n))));
end
n = linspace(0, 101, 101);

%plot
%plotting the next 100 values for each a value --> this should be it's final behaviour.
%If period 1 --> will produce the same value each iteration
%If period 2 --> will produce the same two values each iteration
%If chaotic --> will produce different values each iteration
 

%plots a value against x value
plot(a, x, '.'); 
  hold on; 

[Muhammad Sinan]

Dear @Ruth Porter, you just put the end command at the end of the code, and the code will run. I have just corrected it. So, here is the code while I have tried to run the code here but it took more time than 1 minute therefore it won't run.

%Final plot
%Parameters -- Keep as constants for now.
d = 0.005;
%Parameters - Keep values of b and r constant.
b=8/3;
r=28;
for a = 0.2:0.001:1
%Initial conditions
xold=1.2;
yold=1.3;
zold=1.6;
for n=1:10000000
    xnew=xold+((d*a)*(xold-yold));
    ynew=yold+(d*((r*xold)-(zold*xold)-yold));
    znew=zold+(d*((xold*yold)-(b*zold)));
    xold=xnew;
    yold=ynew;
    zold=znew;
end
x(1)=xnew;
y(1)=ynew;
z(1)=znew;
for n=1:100
    x(n+1)=x(n)+(d*a)*(x(n)-y(n));
    y(n+1)=y(n)+(d)*((r*x(n))-(z(n)*x(n))-y(n));
    z(n+1)=z(n)+(d*((x(n)*y(n))-(b*z(n))));
end
n = linspace(0, 101, 101);
plot(a, x, '.'); 
  hold on; 
end