function dY = SIP(t,init,p)

%% State Variables
Sr = init(1); %initial susceptible red
Sg = init(2); %initial susceptible green
Ir = init(3); %initial infected red
Ig = init(4); %initial infected green
P = init(5); %initial RLO

%% Parameters
Rr = p(1); %
Rg = p(2); %
fr = p(3); %
fg= p(4); %
ur = p(5); %
ug = p(6); %
c = p(7); %
beta = p(8); %
alpha.r = p(9); %
alpha.g = p(10); %
sigma.r = p(11); %
sigma.g = p(12); %
delta = p(13); %

%% ODE
dSr = Rr + fr*(Sr+Ir) - (ur+c*(Sr+Sg+Ir+Ig))*Sr - beta*Sr*P;
dSg = Rg + fg*(Sg+Ig) - (ug+c*(Sr+Sg+Ir+Ig))*Sg - beta*Sg*P;
dIr = beta*Sr*P - (ur+c*(Sr+Sg+Ir+Ig))*Ir - alpha.r*Ir;
dIg = beta*Sg*P - (ug+c*(Sr+Sg+Ir+Ig))*Ig - alpha.g*Ig;
dP = sigma.r*Ir + sigma.g*Ig - delta*P;

%% Output
dY 	= [dSr;dSg;dIr;dIg;dP];