function dy = system_no_recovered(t,y,g_m1, g_f1, g_x1, g_y1, g_m2, g_f2, g_y2, g_x2, beta_mf, beta_fm, beta_fy, beta_yf, beta_yx, beta_xy, beta_mx, beta_xm, eta_m, eta_f, eta_y, eta_x, eta_im, eta_if, eta_iy, eta_ix, n0, alpha_f, alpha_x, alpha_m, alpha_y, alpha_if, alpha_ix, alpha_im, alpha_iy, c_m, c_f, c_y, c_x,mu)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code implements the HIV and Choice Disability model with linear rates of transfer between two choice classes.
% The parameter values are as follows
% t -- time included for compatibility with ODE solvers.  Not used in computation at this time.
% y -- model state.  Should be an 8x1 vector with components in the following order
		% enabled men (m) HIV-
		% enabled women (f) HIV-
		% disabled men (y) HIV-
		% disabled women (x) HIV-
		% enabled men (mi) HIV+
		% enabled women (fi) HIV+
		% disabled men (yi) HIV+
		% disabled women (xi) HIV+
% g -- encode the activity patterns.  Must have g_j1+g_j2=1 for all j=m,f,y,x
% beta -- transmission probabilities per contact.
% eta -- fraction of the incoming population that enters each compartment
% n0 -- total yearly incoming population 
% alpha -- rates of transfer between choice classes. 
		% alpha_m and alpha_f give rate of declining status
		% alpha_x and alpha_y give increasing status
			% similar for infected compartments
% c -- yearly number of contacts for each sex and choice class.
% mu -- univarsal leaving rate.  The current code does not allow for different leaving rates in different choice or infection classes.
% Code originally written by Jeff Musgrave, updated by Rebecca deBoer on 2016/01/26
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			
% death rates

mu_m = mu; mu_f = mu; mu_x = mu; mu_y = mu; mu_im = mu; mu_if = mu; mu_ix = mu; mu_iy = mu; 

% fraction of contacts in each group  According to proportionate mixing in each activity level.

dy = zeros(8,1);

p_1fm = c_m*g_m1*(y(1) + y(5))/(c_m*g_m1*(y(1) + y(5)) + c_f*g_f1*(y(2) + y(6)) + c_y*g_y1*(y(3) + y(7)) + c_x*g_x1*(y(4) + y(8))); %For fixed j, p(1)_ij is the same for all groups

p_1mf = c_f*g_f1*(y(2) + y(6))/(c_m*g_m1*(y(1) + y(5)) + c_f*g_f1*(y(2) + y(6)) + c_y*g_y1*(y(3) + y(7)) + c_x*g_x1*(y(4) + y(8)));

p_1fy = c_y*g_y1*(y(3) + y(7))/(c_m*g_m1*(y(1) + y(5)) + c_f*g_f1*(y(2) + y(6)) + c_y*g_y1*(y(3) + y(7)) + c_x*g_x1*(y(4) + y(8)));

p_1mx = c_x*g_x1*(y(4) + y(8))/(c_m*g_m1*(y(1) + y(5)) + c_f*g_f1*(y(2) + y(6)) + c_y*g_y1*(y(3) + y(7)) + c_x*g_x1*(y(4) + y(8)));

p_1yf = p_1mf; p_1yx = p_1mx; p_1xm = p_1fm; p_1xy = p_1fy;

p_2fm = c_m*g_m2*(y(1) + y(5))/(c_x*g_x2*(y(4) + y(8)) + c_f*g_f2*(y(2) + y(6)) + c_y*g_y2*(y(3) + y(7)) + c_m*g_m2*(y(1) + y(5))); 

p_2yf = c_f*g_f2*(y(2) + y(6))/(c_x*g_x2*(y(4) + y(8)) + c_f*g_f2*(y(2) + y(6)) + c_y*g_y2*(y(3) + y(7)) + c_m*g_m2*(y(1) + y(5))); 

p_2fy = c_y*g_y2*(y(3) + y(7))/(c_x*g_x2*(y(4) + y(8)) + c_f*g_f2*(y(2) + y(6)) + c_y*g_y2*(y(3) + y(7)) + c_m*g_m2*(y(1) + y(5))); 

p_2yx = c_x*g_x2*(y(4) + y(8))/(c_x*g_x2*(y(4) + y(8)) + c_f*g_f2*(y(2) + y(6)) + c_y*g_y2*(y(3) + y(7)) + c_m*g_m2*(y(1) + y(5))); 

p_2xm = p_2fm; p_2mf = p_2yf; p_2xy = p_2fy; p_2mx = p_2yx; 


%%% check for symmetry of contact matrix

%         C_mf = c_m*(y(1) + y(5))*g_m1*p_1mf
%         C_fm = c_f*(y(2) + y(6))*g_f1*p_1fm
%         
%         C_mx = c_m*(y(1) + y(5))*g_m1*p_1mx
%         C_xm = c_x*(y(4) + y(8))*g_x1*p_1xm
%         
%         C_fy = c_f*(y(2) + y(6))*g_f1*p_1fy
%         C_yf = c_y*(y(3) + y(7))*g_y1*p_1yf
%         
%         C_yx = c_y*(y(3) + y(7))*g_y1*p_1yx
%         C_xy = c_x*(y(4) + y(8))*g_x1*p_1xy
%         
%         C_mf = c_m*(y(1) + y(5))*g_m2*p_2mf
%         C_fm = c_f*(y(2) + y(6))*g_f2*p_2fm
%         
%         C_mx = c_m*(y(1) + y(5))*g_m2*p_2mx
%         C_xm = c_x*(y(4) + y(8))*g_x2*p_2xm
%         
%         C_fy = c_f*(y(2) + y(6))*g_f2*p_2fy
%         C_yf = c_y*(y(3) + y(7))*g_y2*p_2yf
%         
%         C_yx = c_y*(y(3) + y(7))*g_y2*p_2yx
%         C_xy = c_x*(y(4) + y(8))*g_x2*p_2xy       

inf=zeros(4,1);
inf(1)=c_m*y(1)*( g_m1*( p_1mf*beta_mf*y(6)/(y(2) + y(6)) + p_1mx*beta_mx*y(8)/(y(4) + y(8))) + g_m2*( p_2mf*beta_mf*y(6)/(y(2) + y(6)) + p_2mx*beta_mx*y(8)/(y(4) + y(8)) ) );
inf(2)=c_f*y(2)*( g_f1*( p_1fm*beta_fm*y(5)/(y(1) + y(5)) + p_1fy*beta_fy*y(7)/(y(3) + y(7)) ) + g_f2*( p_2fm*beta_fm*y(5)/(y(1) + y(5)) + p_2fy*beta_fy*y(7)/(y(3) + y(7)) ) );
inf(3)=c_y*y(3)*( g_y1*( p_1yf*beta_yf*y(6)/(y(2) + y(6)) + p_1yx*beta_yx*y(8)/(y(4) + y(8)) ) + g_y2*( p_2yf*beta_yf*y(6)/(y(2) + y(6)) + p_2yx*beta_yx*y(8)/(y(4) + y(8)) ) );
inf(4)=c_x*y(4)*( g_x1*( p_1xm*beta_xm*y(5)/(y(1) + y(5)) + p_1xy*beta_xy*y(7)/(y(3) + y(7)) ) + g_x2*( p_2xm*beta_xm*y(5)/(y(1) + y(5)) + p_2xy*beta_xy*y(7)/(y(3) + y(7)) ) );

fallm=alpha_m*y(1);
fallf=alpha_f*y(2);
risey=alpha_y*y(3);
risex=alpha_x*y(4);

fallmi=alpha_im*y(5);
fallfi=alpha_if*y(6);
riseyi=alpha_iy*y(7);
risexi=alpha_ix*y(8);

dy(1) = eta_m*n0 - inf(1) - (mu_m)*y(1) - fallm + risey;

dy(2) = eta_f*n0 - inf(2) - (mu_f)*y(2) - fallf + risex;

dy(3) = eta_y*n0 - inf(3) - (mu_y)*y(3) - risey + fallm;

dy(4) = eta_x*n0 - inf(4) - (mu_x)*y(4) - risex + fallf;

dy(5) = eta_im*n0 +inf(1) - (mu_im)*y(5) - fallmi + riseyi;

dy(6) = eta_if*n0 +inf(2) - (mu_if)*y(6) - fallfi + risexi;

dy(7) = eta_iy*n0 +inf(3) - (mu_iy)*y(7) - riseyi + fallmi;

dy(8) = eta_ix*n0 +inf(4)  - (mu_ix)*y(8) - risexi + fallfi;


%c_y*y(3)*( g_y2*( p_2yf*beta_yf*y(6)/(y(2) + y(6))  ) )

%c_y*y(3)*p_2yx*beta_yx*y(8)/(y(4) + y(8))

%c_y*y(3)*( g_y1*( p_1yf*beta_yf*y(6)/(y(2) + y(6)) + p_1yx*beta_yx*y(8)/(y(4) + y(8)) ) ) 

% c_x*y(4)*( g_x1*( p_1xm*beta_xm*y(5)/(y(1) + y(5)) + p_1xy*beta_xy*y(7)/(y(3) + y(7)) ) )
% 
% c_x*y(4)*(g_x2*( p_2xm*beta_xm*y(5)/(y(1) + y(5)) + p_2xy*beta_xy*y(7)/(y(3) + y(7)) ) )

%c_x*y(4)*( g_x2*( p_2xm*beta_xm*y(5)/(y(1) + y(5))  ) )

%a_1 = c_x*y(4)*( g_x1*( p_1xm*beta_xm*y(5)/(y(1) + y(5)) + p_1xy*beta_xy*y(7)/(y(3) + y(7)) ) )

%a_2 = c_x*y(4)*( g_x2*( p_2xm*beta_xm*y(5)/(y(1) + y(5)) + p_2xy*beta_xy*y(7)/(y(3) + y(7)) ) ) 

%c_m*y(1)*( g_m1*( p_1mf*beta_mf*y(6)/(y(2) + y(6)) + p_1mx*beta_mx*y(8)/(y(4) + y(8))) )

%c_m*y(1)*(g_m2*( p_2mf*beta_mf*y(6)/(y(2) + y(6)) + p_2mx*beta_mx*y(8)/(y(4) + y(8)) ) )

%c_y*y(3)*( g_y2*( p_2yf*beta_yf*y(6)/(y(2) + y(6)) + p_2yx*beta_yx*y(8)/(y(4) + y(8)) ) )

%dy(9) = n0 - mu_m*y(9);

%c_m*(y(1) + y(5)) + c_f*(y(2) + y(6)) + c_y*(y(3) + y(7)) + c_x*(y(4) + y(8));