function aids2()

% Predicts epidemic of AIDS
% W is an intermediate individual who has been infected by AIDS but is not
% yet able to transmit disease: represents time delay
% Y is an individual who has been infected and able to transmit disease
% X is susceptibles
% subscript 1 indicates sexually inactive and 2 indicates active population
% F stands for female population of W,Y, and X and M stands for male
% population
% k = average number of sexual contacts in a month
% a = probabilityof contacting AIDS through a single contact
% lamda = immigration rate
% mu = migration rate
% gamma = death rate of AIDS
% p = high contact fraction
% tau = ratio of actives to inactives
% b = rate of time delay from W to Y
% q = ratio of female to male

% Total numbers of W,Y and X are calculated by adding F and M

% For female:
% dWf1/dt = kaf*Xf1*(Ym1+tau*Ym2)/(Xf1+Yf1+tau*(Yf2+Xf2)) - (b+mu)*Wf1
% dWf2/dt = kaf*tau*Xf2*(Ym1+tau*Ym2)/(Xf1+Yf1+tau*(Yf2+Xf2)) - (b+mu)*Wf2
% dYf1/dt = b*Wf1 - (mu+gamma)*Yf1
% dYf2/dt = b*Wf2 - (mu+gamma)*Yf2
% dXf1/dt = -kaf*Xf1*(Ym1+tau*Ym2)/(Xf1+Yf1+tau*(Yf2+Xf2)) + (1-p)*.5*lamda - mu*Xf1
% dXf2/dt = -kaf*tau*Xf2*(Ym1+tau*Ym2)/(Xf1+Yf1+tau*(Yf2+Xf2)) + p*.5*lamda - mu*Xf2
% For male:
% dWm1/dt = ka*Xm1*(Ym1+tau*Ym2)/(Xm1+Ym1+tau*(Ym2+Xm2)) - (b+mu)*Wm1
% dWm2/dt = ka*tau*Xm2*(Ym1+tau*Ym2)/(Xm1+Ym1+tau*(Ym2+Xm2)) - (b+mu)*Wm2
% dYm1/dt = b*Wm1 - (mu+gamma)*Ym1
% dYm2/dt = b*Wm2 - (mu+gamma)*Ym2
% dXm1/dt = -ka*Xm1*(Ym1+tau*Ym2)/(Xm1+Ym1+tau*(Ym2+Xm2)) + (1-p)*.5*lamda - mu*Xm1
% dXm2/dt = -ka*tau*Xm2*(Ym1+tau*Ym2)/(Xm1+Ym1+tau*(Ym2+Xm2)) + p*.5*lamda - mu*Xm2

% time variables
tspan = [0:1:12*5];

% initial conditions
% female infective is initially 18% of infective population
% data from 1997 is used
% it was assumed that 90 % of infective population is sexually active
% [Wf1, Wf2, Yf1, Yf2, Xf1, Xf2, Wm1, Wm2, Ym1, Ym2, Xm1, Xm2]
Wm1=5000; Wm2=40000; Ym1=round(641086*.82*.1-5000); Ym2=round(641086*.82*.9-40000); Xm1=102900000; Xm2=44100000;
Wf1=round(.1*Wm1); Wf2=round(.1*Wm2); Yf1=round(641086*.18*.1-Wf1); Yf2=round(641086*.18*.9-Wf2); Xf1=Xm1; Xf2=Xm2;
D = 0;      % Calculates cumulative death
init = [Wf1, Wf2, Yf1, Yf2, Xf1, Xf2, Wm1, Wm2, Ym1, Ym2, Xm1, Xm2, D];

% parameters
ka = .023;
b = 1/(12*2);
mu = .00417;
lamda  = mu*(Xf1+Xf2+Xm1+Xm2);
tau = 5;
gamma = 1/(12*3);
p = .2;
kaf = .001;

% resolution of solution can be changed by changing refine factor
options = odeset('Refine', 1);

%commented for loops can be used to generate a plot with varying parameter
%values
%for i = 1:4
[t,u] = ode45(@ODE, tspan, init, options, ka, b, mu, lamda, tau, gamma, p, kaf);
%S3 = (u(:,3)+u(:,4)+u(:,9)+u(:,10))./(u(:,1)+u(:,2)+u(:,3)+u(:,4)+u(:,5)+u(:,6)+u(:,7)+u(:,8)+u(:,9)+u(:,10)+u(:,11)+u(:,12));
%ka = ka - .003*i;
%A(:,i) = S3;
%end

%Total number of AIDS infective at the end of tspan
T = u(length(u),3)+u(length(u),4)+u(length(u),9)+u(length(u),10)

% codes used to plot the result of for loop
%plot(t, A(:,1), 'm', t, A(:,2), 'b',t, A(:,3), 'y', t, A(:,4),'g');
%legend('ka=.023', 'ka=.020', 'ka=.014', 'ka=.005')
%xlabel('time (months)'); ylabel('Fraction of Infective'); title('Fraction of Infective with varying ka values');
%figure;

%plot cumulative death
%plot(t, u(:,13));

%calculates fraction of infective for female, male and total respectively
%S1 = (u(:,3)+u(:,4))./(u(:,1)+u(:,2)+u(:,3)+u(:,4)+u(:,5)+u(:,6));
%S2 = (u(:,9)+u(:,10))./(u(:,7)+u(:,8)+u(:,9)+u(:,10)+u(:,11)+u(:,12));
%S3 = (u(:,3)+u(:,4)+u(:,9)+u(:,10))./(u(:,1)+u(:,2)+u(:,3)+u(:,4)+u(:,5)+u(:,6)+u(:,7)+u(:,8)+u(:,9)+u(:,10)+u(:,11)+u(:,12));

% plots intermediates and infectives
figure;
plot(t, u(:,1)+u(:,7),'bo');hold on;
plot(t, u(:,2)+u(:,8),'bx');hold on;
plot(t, u(:,3)+u(:,9),'g+');hold on;
plot(t, u(:,4)+u(:,10),'gs');
legend('W1','W2','Y1','Y2');xlabel('time (months)');ylabel('Population');title('Number of Intermediates and Infectives');

% plots susceptibles
figure;
plot(t, u(:,5)+u(:,11),'md');hold on;
plot(t, u(:,6)+u(:,12),'mv');
legend('X1','X2'); xlabel('time (months)'); ylabel('Population');title('Number of Susceptibles');




function uprime = ODE(t,u,ka,b,mu,lamda,tau,gamma,p,kaf)
% ODEs
% u(1) = Wf1, u(2) = Wf2, u(3) = Yf1, u(4) = Yf2, u(5) = Xf1, u(6) = Xf2
% u(7) = Wm1, u(8) = Wm2, u(9) = Ym1, u(10) = Ym2, u(11) = Xm1, u(12) = Xm2

uprime = [((ka2*u(5)*(u(9)+tau*u(10)))/(u(5)+u(3)+tau*(u(4)+u(6))) - (b+mu)*u(1));
          ((ka2*tau*u(6)*(u(9)+tau*u(10)))/(u(5)+u(3)+tau*(u(4)+u(6))) - (b+mu)*u(2));
          (b*u(1) - (mu+gamma)*u(3));
          (b*u(2) - (mu+gamma)*u(4));
          (-(ka2*u(5)*(u(9)+tau*u(10)))/(u(5)+u(3)+tau*(u(4)+u(6))) + (1-p)*.5*lamda - mu*u(5));
          (-(ka2*tau*u(6)*(u(9)+tau*u(10)))/(u(5)+u(3)+tau*(u(4)+u(6))) + p*.5*lamda - mu*u(6));
          ((ka*u(11)*(u(9)+tau*u(10)))/(u(11)+u(9)+tau*(u(10)+u(12))) - (b+mu)*u(7));
          ((ka*tau*u(12)*(u(9)+tau*u(10)))/(u(11)+u(9)+tau*(u(10)+u(12))) - (b+mu)*u(8));
          (b*u(7) - (mu+gamma)*u(9));
          (b*u(8) - (mu+gamma)*u(10));
          (-(ka*u(11)*(u(9)+tau*u(10)))/(u(11)+u(9)+tau*(u(10)+u(12))) + (1-p)*.5*lamda - mu*u(11));
          (-(ka*tau*u(12)*(u(9)+tau*u(10)))/(u(11)+u(9)+tau*(u(10)+u(12))) + p*.5*lamda - mu*u(12))
          gamma*(u(3)+u(4)+u(9)+u(10))];