clear all
x(1)=0.3;  %the initial value of the sequence
N=500; % number of iterations
lambda=3.3; %lambda is the maximum growth rate per capita
beta=4; % parameter beta
z=[x(1) x(1)]; % stores the pairs of points to plot
for n=1:N-1  % begins the loop
    x(n+1)=lambda*x(n)/(1+x(n)^beta); % iterate the difference equation
    z=[z; x(n) x(n+1); x(n+1) x(n+1)];  % stores points
end

% the following lines format the plot in a nice way
a=min(x)-(max(x)-min(x))/10;
a=0; b=max(x)+(max(x)-min(x))/10;
delx=(b-a)/100; x=[a:delx:b]; 
plot(x,lambda.*x./(1+x.^beta),'-r',x,x,'-k','LineWidth',2)
hold on
plot(z(:,1),z(:,2),'-b','LineWidth',1)
xlim([a,b])
ylim([a,b])
hold off
