function [pCount_X pCount_Y pCount_Z] = discreteSIR(N,beta,gamma,T,K)
% discrete-time, discrete-space, run 200times
% if a individual is a susceptible, the value is 0, infected 1, recovery 2;
rand('state',sum(100*clock)); %random seed
%Initialize
for k = 1 : K
Population_List = zeros(1,N);
P_Zero = randperm(100);%random 100 infected individuals
for i = 1 : size(P_Zero,2)
Population_List(P_Zero(i)) = 1;
end
Count_X{k} = []; % calculate the susceptible individuals in each time step;
Count_Y{k} = []; % calculate the infected individuals in each time step;
Count_Z{k} = []; % calculate the recovery individuals in each time step;
pCount_X{k} = []; %probability of X;
pCount_Y{k} = []; %probability of Y;
pCount_Z{k} = []; %probability of Z;
new_Y{k} = [];%the newly cases in each time step
new_Y{k}(1,1) = 0;
new_R{k} = [];%the newly recovery in each time step
new_R{k}(1,1) = 0;
Count_Y{k}(1,1) = sum(Population_List == 1);
Count_X{k}(1,1) = N - Count_Y{k}(1,1);
while sum(Population_List == 1) > 0
sumX = 0;
sumY = 0;
temp_X = [];
temp_X = find(Population_List == 0);
for i = 1 : size(temp_X,2)
if rand < 0.4 && (Population_List(1 + floor(rand() * N)) == 1)
Population_List(temp_X(i)) = 1;
sumX = sumX + 1;
end
end
new_Y{k} = [new_Y{k}, sumX];
temp_Y = [];
temp_Y = find(Population_List == 1);
for j = 1 : size(temp_Y,2)
if rand < gamma
Population_List(temp_Y(j)) = 2;
sumY = sumY + 1;
end
end
new_R{k} = [new_R{k}, sumY];
end
%calculate the number of Susceptible, Infected and recovery
for i = 2 : T
Count_X{k}(i) = Count_X{k}(i-1) - new_Y{k}(i);
Count_Y{k}(i) = Count_Y{k}(i-1) + new_Y{k}(i)- new_R{k}(i);
end
Count_Z{k} = N - Count_X{k} - Count_Y{k};
pCount_X{k} = Count_X{k} / N;
pCount_Y{k} = Count_Y{k} / N;
pCount_Z{k} = Count_Z{k} / N;
end
%figure
%for k = 1 : 200
% plot(pCount_X{k},'ok','MarkerSize',3);
% hold on
% plot(pCount_Y{k},'ok','MarkerSize',3);
%hold on
%plot(pCount_Z{k},'ok','MarkerSize',3);
%end