%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      Harvesting model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     X(n+1)= lambda X(n) - h
%
% with initial condition X(0)=X0
% has exact solution:
%
%   X(n)=lambda^n X0 -(lambda^n-1)/(lambda-1) h
%
% First calculate X(n) without the formula above.
%
clear X;
lambda=1.1;
X(1)=2;
%h=0.20;
h=0.30;
N=10;
for n=1:N
    X(n+1)=lambda*X(n)-h;
end;
figure(4)
plot(0:N,X)
%
% Now compute Y=X using formula above
%
m=0:N;
Y=lambda.^m*X(1) -(lambda.^m-1)./(lambda-1)*h;
hold on
plot(m,Y,'r*')
title(['x_n in harvest model: lambda=',num2str(lambda),'h=',num2str(h)])
xlabel('n')
ylabel('x_n')
hold off

%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      Harvesting model with Seasonal Variation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     Seasonal variation of growth rate. Suppose 
%     n=time (in months) since start growth rate at 
%     time n is L(n):
%
%     L = L0+alpha*sin(omega*n)
%
%     X(n+1)= L(n) X(n) - h
%
%     with initial condition X(1)=X0.
% 
%
clear X;
N=120;
L0=1.07;
alpha=0.5;
omega=2*pi./12;
L=L0+alpha*sin(omega*(1:N));
%
%
X(1)=2;
h=0.000;
%h=0.035;
for n=1:N
    X(n+1)=L(n)*X(n)-h;
end;
%
figure(5)
clf
hold on
plot(0:N,X,'r-')
plot(1:N,L,'b-')
title(['x_n in seasonal model: L0=',num2str(L0),'  h=',num2str(h)])
xlabel('n')
ylabel('x_n')
axis([0,N,0,20])
hold off

