function P2()
close all;
t = linspace(0, 1, 100);

figure(1);
hold;
Y0 = [95; 5; 0];
c = 1;
d = 5;
[y1 y2 y3] = Kermack_McKendrick(Y0, t, c, d);
plot(t, y1, t, y2, t, y3);
legend('y1', 'y2', 'y3');
title(sprintf('c=%.2f, d=%.2f, y(0)=(%.2f, %.2f, %.2f)', c, d, Y0(1), Y0(2), Y0(3)));

figure(2);
%   y2 does not increase
% => c*y1*y2-d*y2 <= 0
% => c*y1-d <= 0
% Because y1' <= 0, if c*y1(0)-d <= 0, then c*y1(t)-d<=0 for all t>0.
% we let c = 0.1, d = 10, y1(0) = 75, y2(0) = 25, y3(0) = 0.
Y0 = [75; 25; 0];
c = 0.1;
d = 10;
[y1 y2 y3] = Kermack_McKendrick(Y0, t, c, d);
plot(t, y1, t, y2, t, y3);
legend('y1', 'y2', 'y3');
title(sprintf('c=%.2f, d=%.2f, y(0)=(%.2f, %.2f, %.2f)', c, d, Y0(1), Y0(2), Y0(3)));

figure(3);
% we let c = 0.23, d = 5, y1(0) = 95, y2(0) = 5, y3(0) = 0.
Y0 =  [95;5;0];
c = 0.23;
d = 5;
[y1 y2 y3] = Kermack_McKendrick(Y0, t, c, d);
plot(t, y1, t, y2, t, y3);
legend('y1', 'y2', 'y3');
title(sprintf('c=%.2f, d=%.2f, y(0)=(%.2f, %.2f, %.2f)', c, d, Y0(1), Y0(2), Y0(3)));


function [y1 y2 y3] = Kermack_McKendrick(Y0, t, c, d)
% Y0 =  [y10; y20; y30]: initial value
% [y1 y2 y3] are value functions at time t
% ODE: y1' = -c*y1*y2
%      y2' = c*y1*y2-d*y2
%      y3' = d*y2
if size(Y0, 1) == 1
    Y0 = Y0'; % Y0 should be column vector
end
L = length(t);
y = zeros(3,L);
y(:,1) = Y0;
options = odeset;
for i = 2:L
    [time ytemp] = ode23(@f, [t(i-1) t(i)], y(:,i-1),options,[c;d]);
    y(:, i) = ytemp(end,:)';
end
y1 = y(1,:);
y2 = y(2,:);
y3 = y(3,:);

function  dydt = f(t, y, p)
% y: [y1 y2 y3]
% p: [c d]
dydt = zeros(3,1);
dydt(1) = -p(1)*y(1)*y(2);
dydt(2) = p(1)*y(1)*y(2)-p(2)*y(2);
dydt(3) = p(2)*y(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3 points for correct  f function and figure 1 plot
% 1 point for correct figure 2
% 1 point for correct figure 3   something like 1.01  good enough accuracy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%