%cobweb.m
%This is a script to iterate a first order nonlinear map
%You should only need to edit the first block of code


%This controls how many iterations you compute
N = 1000;
%This is where we define parameters
a = 5.9;
b = 2.03;
%this is the inital value of the model
x0 = 2.1;

%This is where we define the map which acts on x
%Comment out one of the following lines (put a % in the front) and 
%the computer will ignore that case

mapping = @(x) x - x.*(x.^a)./(2^a + x.^a) + 1; %Problem 1
% mapping = @(x) a.*x.^(1-b); %Problem 2





%getting ready to run simulations
y = [x0;x0];
z = linspace(0,2.5,100);
fz = mapping(z);
x = x0;
xz = [x];

%We now run forward in time until 'N'
for j = 1:N
    
    %Apply the map to our variable
    f = mapping(x);
    %Get ready for next step
    x = f;
    %add the value to 'xz' for plot #2
    xz = [xz;x];
    %Add the value to 'y' for plot #1
    y = [y;x;x];
    
end

%Here we plot the function, the line x=y, and our cobweb
figure(1)
plot(y(1:2*N+1),[0;y(3:2*N+2)],z,fz,'r',z,z,'b--', 'LineWidth',2)
title('Cobweb plot of system','FontSize',18)
xlabel('Input','Fontsize',16)
ylabel('Output','Fontsize',16)

%Here we will plot the time sequence of our solution
figure(2)
plot(xz,'*')
title('Solution sequence','FontSize',18)
xlabel('n','Fontsize',16)
ylabel('Output','Fontsize',16)