function N=cobweb(A,N0,T,maxN)
% function N=cobweb(A,N0,T,maxN)
%
% A is a parameter vector to be passed on to seqfnc.
%
% This function calculates the sequence up to element T
%  for the difference equation N_(t+1)=seqfnc(N_t) with N_0=N0.
% It plots the cobweb diagram and the sequence N.
%
% maxN is the largest value of N for the plots.


% define the parameters

R = A(1);
K = A(2);

% compute the points to plot on the graph y=f(N)

x = [maxN/100:maxN/100:maxN];
f = 0;
for j=1:length(x)
    f(j)=seqfnc(x(j));
end

% set up the vector of times

t=[0:T];

% compute the sequence values by iteration
% "N=0" is used to erase any values previously stored as "N".
% (xx,yy) are the points (N_t,N_{t+1})
% (webx,weby) are coordinates for the stairstep cobweb line

N = 0;
webx = 0;
weby = 0;
xx = 0;
yy = 0;
N(1)=N0;
webx(1)=N0;
weby(1)=0;
for i=1:T
    N(i+1)=seqfnc(N(i));
    xx(i)=N(i);
    yy(i)=N(i+1);
    webx(2*i)=N(i);
    weby(2*i)=N(i+1);
    if i<T
        webx(2*i+1)=N(i+1);
        weby(2*i+1)=N(i+1);
    end
end

% subplot(1,2,k) indicates plot k, given 1 row and 2 columns

% Plot y=f(N) as a thick blue curve.
% Plot y=N in black and the cobweb segments in red.
% Mark the points (N_t,N_{t+1}) as black diamonds.
% Set the axis ranges and labels.

subplot(1,2,1)
plot(x,f,'b','LineWidth',2)
hold on
plot([0:maxN:maxN],[0:maxN:maxN],'k','LineWidth',1)
plot(webx,weby,'r')
plot(xx,yy,'kd')
xlim([0,maxN])
ylim([0,maxN])
xlabel('\it{N}','FontSize',12)
ylabel('\it{N}','FontSize',12,'Rotation',0)
hold off

% Plot the points (t,N) as "x".  Set the axis ranges and labels.

subplot(1,2,2)
plot(t,N,'kx')
xlim([0,T])
ylim([0,maxN])
xlabel('\it{t}','FontSize',12)
ylabel('\it{N}','FontSize',12,'Rotation',0)

% The code that follows is used to compute the new population levels
% The parameters need to have been defined earlier

    function y=seqfnc(N)
        y = N+R*N*(1-N/K);
    end

end