function seir = seir(t,x)
N = 1000;
beta = 0.25;
gamma = 0.1;
epsilon = 0.2;
seir = zeros(3,1);
seir(1) = - beta * (x(1) / N) * x(3); %dS/dt
seir(2) = beta * (x(1) / N) * x(3) - epsilon * x(2); %dE/dt
seir(3) = epsilon * x(2) - gamma * x(3); %dI/dt

