%Set up passive regulation function
at=[0:0.01:800];
at1=alpha+at.*exp(-at/n);

n=1000;
alpha=100;


a(1)=1;
t=[1:1:200];

%Iterate the regulatory function
for i=t
    a(i+1)=alpha+a(i).*exp(-a(i)/n);
end

%Make a cobweb diagram
j=0;
for i=1:1:(length(a)-1)
    j=j+1;
    a2(j)=a(i);
    a3(j)=a(i);
    j=j+1;
    a2(j)=a(i);
    a3(j)=a(i+1);
end

%Plot cobweb diagram
figure(1)
plot(at,at1,'k-',at,at,'k--',a2,a3,'k:')
xlabel('Population at time t: x_t')
ylabel('Population at time t+1:x_{t+1}')
axis([0 500 0 500])
title('Passive feedback')

%Plot time series
figure(2)
plot(t,a(1:max(t)),'k')
xlabel('Time: t')
ylabel('Population :x_{t}')
axis([0 20 0 500])
title('Passive feedback')



%Set up active regulation function
at=[0:0.01:7500];

n=1000;
b=20;
at1=b*at.*exp(-at/n);

%Iterate the regulatory function
a(1)=1;
t=[1:1:200];
for i=t
    a(i+1)=b*a(i).*exp(-a(i)/n);
end


%Make a cobweb diagram
j=0;
for i=1:1:(length(a)-1)
    j=j+1;
    a2(j)=a(i);
    a3(j)=a(i);
    j=j+1;
    a2(j)=a(i);
    a3(j)=a(i+1);
end

%Plot cobweb diagram
figure(3)
plot(at,at1,'k-',at,at,'k--',a2,a3,'k:')
xlabel('Population at time t: x_t')
ylabel('Population at time t+1:x_{t+1}')
axis([0 7500 0 7500])
title('Active feedback')

%Plot time series
figure(4)
plot(t,a(1:max(t)),'k')
xlabel('Time: t')
ylabel('Population :x_{t}')
axis([0 20 0 7500])
title('Active feedback')