options = odeset('RelTol',1e-4,'AbsTol',[1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9 1e-9]);

N=2000000;
[T,Y] = ode45(@CDrun,[0 70],[0.9 0.0 2.7700e-04 0.0 0.1 0.01 0.9 0.0 2.7700e-04 0.1 0.01],options);


plot(T,Y(:,3)*N, '-',  T,Y(:,9)*N, 'r--', 'MarkerSize',8, 'LineWidth',2.5,'Color','black');
legend('With multiple control measures',' Without control measures','Location','Best');
xlabel('Time (days)', 'FontSize', 14);
ylabel('Number of infectious individuals', 'FontSize', 14);

function dy = CDrun(t,y)
dy = zeros(11,1);    % a column vector

mu = 0.002; % Chen et al. Infectious Diseases of Poverty 
delta  = 0.113; %.High estimated per day
gamma  = 0.172; % Chen et al. Infectious Diseases of Poverty
xi = 0.100;   % Chen et al. Infectious Diseases of Poverty

% parameters fit to data
beta = 1.097; %4828;
alpha =  0.098; %1665;
phi = 0.006; %0096;    
sigma = 0.947;  %9056;
rho = 0.242;  %;
varepsilon = 0.211;

R0c = (mu*sigma*(alpha+(1-varepsilon)*beta))/((mu+phi)*(mu+sigma)*(mu+rho+delta+gamma));

R0 = (sigma*(alpha+beta))/((mu+sigma)*(mu+delta+gamma));


dy(1) = mu-(1-varepsilon)*beta*y(1)*y(3)-alpha*y(1)*y(6)-(mu+phi)*y(1);
dy(2) = (1-varepsilon)*beta*y(1)*y(3)+alpha*y(1)*y(6)-(mu+sigma)*y(2);
dy(3) = sigma*y(2)-(mu +rho + delta + gamma)*y(3);
dy(4) = rho*y(3)-(mu+delta+gamma)*y(4); 
dy(5) = phi*y(1)+(delta+gamma)*(y(3)+y(4))-mu*y(5);
dy(6) = xi*(y(3)-y(6));

dy(7) = mu-beta*y(7)*y(9)-alpha*y(7)*y(11)-mu*y(7);
dy(8) = beta*y(7)*y(9)+alpha*y(7)*y(11)-(mu+sigma)*y(8);
dy(9) = sigma*y(8)-(mu + delta + gamma)*y(9);
dy(10) =(delta+gamma)*y(9)-mu*y(10);
dy(11) = xi*(y(9)-y(11));
end