clear all;close all;
Nx = 1000; 
Ntherm = 1000;
rt = 2.999:0.00001:4;
res=[];
rbif=[];
cold = 1;
for ir = 1:size(rt,2)
    r = rt(ir);
    %disp(['r=' num2str(r)]);
    x = 0.6323424;
    xt=[];
    c = 0;
    for iter=-Ntherm:Nx
        x = r.*x.*(1-x);
        if iter>0
            if c==0
                for j = 1: iter-1
                    if abs(x-xt(j))<10^-6;c=iter-1;end
                end
            end
            xt = [xt x];
        end 
    end
    if c>cold 
        rbif = [rbif r]
    end
    cold = c;
    xt = sort(xt);
    res = [res [r xt]'] ;     
    %plot(res(1,ir).*ones(Nx,1),res(2:Nx+1,ir),'b.','MarkerSize',1);hold on; 
    rbif;
   
end

for k = 2 : size(rbif,2)-1
    feig = (rbif(k)-rbif(k-1))/(rbif(k+1)-rbif(k));
    disp([ 'k=' num2str(k) ' Feigenbaum constant = ' num2str(feig)])
end
