clear('x')
N=50;
r=3.3;
b=4;
x(1)=b/2;
for n=1:N
    x(n+1)=r*x(n)/(1+x(n)^b)
end
plot([1:1:N+1],x,'-o')
