
function dy = func_SIS(t,y)

global alpha beta0 beta1 gamma mu

S = y(1);
E = y(2);
I = y(3);
tt = y(4);

dS = mu - beta_SEIR(tt)*S*I - mu*S;
dE = beta_SEIR(tt)*S*I - (mu + alpha)*E;
dI = alpha*E - (mu+gamma)*I;
ds = 1;

dy = [dS;dE;dI;ds];
