% States X
% X(1) = Uninfected
% X(2) = Infected
% X(3) = Cured
% X(4) = Dead;
X = [700000;100;0;0];

CureRate = 0.05;
DeathRate = 0.007;
Infectivity = .0348;
N = 10;             % interractions per person

Y = [];

for i=1:300
    
    NewInfections = X(1) * X(2) * Infectivity * N / sum(X(1:3)) ;
    Cures = CureRate * X(2);
    Deaths = DeathRate * X(2);
    
    X(1) = X(1) - NewInfections
    X(2) = X(2) + NewInfections - Cures - Deaths;
    X(3) = X(3) + Cures;
    X(4) = X(4) + Deaths;

    Y = [Y ; X'];
end

plot(Y/1000)
xlabel('days');
ylabel('x1000');
title('Healthy(b), Infected(g), Cured(r), Dead(c)');