%
% cobwebbing for Xm+1 = f(Xn)
%

small_k = input('input small k = ? ');
N = input('saturation level = ? ');
r = 1 + small_k*N;
K = r/small_k;

x = linspace(0, K);
fx = r*(1-x/K).*x;
plot(x,fx), pause

hold on
plot(x,x,'-m'),
pause

x0 = input('initial population = ? '); 
for i = 1:1000
    x1 = r*(1-x0/K)*x0;
    z = [x0 x0]; w = [x0 x1];
    plot(z,w,'Or-'), pause
    z = [x0 x1]; w = [x1 x1];
    plot(z,w,'Oc-'),pause
    x0 = x1;
end