function aids()

% Predicts epidemics 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 male to female

% Total number of W,Y and X are calculated by adding F and M
% Assumptions:
% 1) F and M have the same rates of following: k, lamda,mu, and p
% 2) Male population has higher rates of obtaining AIDS.
%    In other words, initial population of susceptibles and infectives of
%    male are greater than those of female

% Both Male and Female have the same equations except the different values
% of rates
% dW1/dt = k*a*X1*(Y1+tau*Y2)/(X1+Y1+tau*(Y2+X2)) - (b+mu)*W1
% dW2/dt = k*a*tau*X2*(Y1+tau*Y2)/(X1+Y1+tau*(Y2+X2)) - (b+mu)*W2
% dY1/dt = b*W1 - (mu+gamma)*Y1
% dY2/dt = b*W2 - (mu+gamma)*Y2
% dX1/dt = -k*a*X1*(Y1+tau*Y2)/(X1+Y1+tau*(Y2+X2)) + (1-p)*lamda - mu*X1
% dX2/dt = -k*a*tau*X2*(Y1+tau*Y2)/(X1+Y1+tau*(Y2+X2)) + p*lamda - mu*X2

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

% initial conditions
% [W1, W2, Y1, Y2, X1, X2]

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;

init = [Wm1+Wf1,Wm2+Wf2,Ym1+Yf1,Ym2+Yf2,Xm1+Xf1,Xm2+Xf2];
% Female population parameters
ka = .05;
b = 1/(12*2);
mu = .00417;
lamda  = mu*(Xf1+Xf2+Xm1+Xm2);
tau = 5;
gamma = 1/(12*3);
p = .2;

options = odeset('Refine', 1);
[t,F] = ode45(@ODE, tspan, init, options, ka, b, mu, lamda, tau, gamma, p);

S = (F(:,3)+F(:,4))./(F(:,1)+F(:,2)+F(:,3)+F(:,4)+F(:,5)+F(:,6));
F(length(F),3)+F(length(F),4)

F(length(F),3)+F(length(F),4);
D = mu.*(F(:,3)+F(:,4));
figure;
plot(t, D)
figure;
plot(t,F(:,2),'bx',t,F(:,3),'g+',t,F(:,4),'gs',t,F(:,5),'md',t,F(:,6),'mv');
figure;
plot(t, S)

function uprime = ODE(t,u,ka,b,mu,lamda,tau,gamma,p)
% ODEs
% u(1) = W1, u(2) = W2, u(3) = Y1, u(4) = Y2, u(5) = X1, u(6) = X2

uprime = [((ka*u(5)*(u(3)+tau*u(4)))/(u(5)+u(3)+tau*(u(4)+u(6))) -(b+mu)*u(1));
          ((ka*tau*u(6)*(u(3)+tau*u(4)))/(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));
          (-(ka*u(5)*(u(3)+tau*u(4)))/(u(5)+u(3)+tau*(u(4)+u(6))) + (1-p)*lamda - mu*u(5));
          (-(ka*tau*u(6)*(u(3)+tau*u(4)))/(u(5)+u(3)+tau*(u(4)+u(6))) + p*lamda - mu*u(6))];