calculate orbit diagram of the logistic map

close;
clear;
orbit=zeros(1,300);
j=0;
% parameter range
for(r=2.8:0.001:4)
    j=j+1;
    % random nitiation of iteration
    xn1=rand(1);
        for(i=1:600)
        % calculate logistic map
        xn=xn1;
        xn1=r*xn*(1-xn);
       % wait for transients
       if(i>300)
           % store the orbit points
           orbit(i-300)=xn1;
       end
    end
    plot(r,orbit);
    if(j==1)
        axis([2.8 4 0 1]);
    hold;
    end
    end
Current plot held