function [ydot] = epidem(t,y)
    % function [ydot] = epidem(t,y)
    % t is scalar for time 
    % y = [S,D,M,I] is vector for unknown functions 
    % ydot is vector of right-hand-sides of the system of differential equations
    % a, b, r are parameters of the system
    
a = 6.90675*10^(-5);
b = 10;
r = 0.1;
ydot = [ -a*y(1)*y(4); r*y(4); b; a*y(1)*y(4) - r*y(4) - b]; 