% Solution to problem 1 of control 1A
% JMS. Oct.2020

% Set parameters
clear
version = 'A'
if version=='A'
  a = 3.0;
  b = 4e-2;
  c = 0.5;
  d = 1e-3;
  x1 = 700;
  y1 = 80;
elseif version=='B'
  a = 1.0;
  b = 1e-2;
  c = 0.8;
  d = 4e-4;
  x1 = 2666;
  y1 = 333;
elseif version=='C'
  a = 1.5;
  b = 4e-3;
  c = 0.4;
  d = 3e-4;
  x1 = 1500;
  y1 = 250;
elseif version=='D'
  a = 4.0;
  b = 8e-3;
  c = 0.7;
  d = 6e-4;
  x1 = 2000;
  y1 = 500;
end
h = 5e-4;
tmax = 20;

% Find populations
A = [a,0;0,-c];
B = [0,-b;d,0];
t = 0:h:tmax;
nt = numel(t);
q = [x1;y1];
for it = 1:nt-1
  q(:,it+1) = (eye(2)+h*A)*q(:,it) +h*(B*q(:,it)).*q(:,it);
end
x = q(1,:);
y = q(2,:);

% Plot populations
figure(11)
plot(t,x, t,y, 'LineWidth',2)
legend('preys','predators')
xlabel('time (yr)')
ylabel('population')
set(gca,'FontSize',20)
grid on

% Print requested times
fprintf('               Time with x>x1 = %4.1f yr\n',sum(x<x1)*h)
fprintf('     Time with x<x1 or x>2*x1 = %4.1f yr\n',sum((x<x1)|(x>2*x1))*h)
fprintf('First time with x<x1 and y<y1 = %4.1f yr\n',find((x<x1)&(y<y1),1)*h)

% Plot predator population with colour according to rate
figure(12)
it = 1:nt-1;
mu = (y(it+1)-y(it))./(h*y(it));
mu0 = 0.01;
ir = mu>+mu0;
ig = mu<-mu0;
ib = -mu0<=mu & mu<=+mu0;
t = t(it);
y = y(it);
plot(t(ir),y(ir),'r.',t(ig),y(ig),'g.',t(ib),y(ib),'b.')
legend('\mu>0.01','\mu<-0.01','-0.01<=\mu<=0.01')
xlabel('time (yr)')
ylabel('predators')
set(gca,'FontSize',20)
grid on
