function r=Feigenbaum(N)
%Finds sequence of superstable orbits and computes delta adn alpha

    crit = 0.5;     %critical point
 
%Find the period 1 and 2
    r(1)=superstable(1,2,crit);
    r(2)=superstable(2,3,crit);
    d(1)=0;
    d(2)= Logistic(1,r(2),crit)-crit;
    delta(1)=4.0;
    delta(2)=4.0;
%Loop to find superstable period 2^(i-1)
    for i=3:N
   
        period =2^(i-1);
        rguess= r(i-1) + (r(i-1)-r(i-2))/delta(i-1);        %Use Geometric scaling for guess        
        rold = 0.7*rguess+0.3*r(i-1);                       %need a second guess for secant

        r(i) = superstable(period,  rguess,rold);           %find orbit
 
        d(i) = Logistic(period/2,r(i), crit) - crit;		%compute width of branch
        %Feigenbaum scaling parameters
        delta(i) = (r(i-1)-r(i-2))/(r(i)-r(i-1));
        alpha(i) = (d(i-1)-d(i-2))/(d(i)-d(i-1));
    end
    
    fprintf('period\t\t  r_ss\t\t  delta\t\talpha\n');
    for i=1:N
         fprintf('%d\t%.14f\t%.11f\t%.11f\n', 2^(i-1), r(i), delta(i),alpha(i));   
    end
    