diff --git a/ODEsImmunity.R b/ODEsImmunity.R
index b901009..14fb4f5 100644
--- a/ODEsImmunity.R
+++ b/ODEsImmunity.R
@@ -23,7 +23,7 @@ SEAI<-function(t,y,p){
     dPSU.dt = -(b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSU)*(PSU)
     dPSD.dt = -(b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSU/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSD)*(PSD)
     dPRU.dt = -(b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRU)*(PRU)
-    dPRD.dt = -(b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRU/N)*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRD)*(PRD)
+    dPRD.dt = -(b*(PRU/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRD)*(PRD)
     #dNU.dt = (del*(AU+AD))+(ND/Tau)-(b*I*eps*NU*(d+rho1))-(muNU*NU) - (del*att*NU)
     dNU.dt = del+(ND/Tau)-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU)
     dND.dt = (b*I*eps*NU*d1)-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ND*I*b*eps) - (del*att*ND)
@@ -37,13 +37,13 @@ SEAI<-function(t,y,p){
   })
 }
 
-b = 0.9
+b = 0.4*0.01
 
 nstrains = 20
 
 att = 0.001#.000000000001
 del = 60
-muND = 0.05
+muND = 0.1
 muNU = 0.1
 rho1 = 0.5 #prob of immunity gain from NU to SU
 rho2 = 0.5 #prob of immunity gain from ND to SD
@@ -58,19 +58,19 @@ L = 10000
 
 
 ws = 1
-wr = 1-0.6#-0.6
+wr = 1-0.6
 #muSU = 1
-muSU=0.05
+muSU=0.001
 #muSD = 20
 muSD=0.99
 #muRD = 1
-muRD = 0.05
+muRD = muSU#0.001
 #muRU = 1
-muRU = 0.05
+muRU = muSU#0.001
 Tau = 20#/365
 d1 = 0.3
-d2 = 0.01
-d3 = 0.01
+d2 = d1/6
+d3 = d1/10
 #m = 3
 K=10
 p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att)
@@ -79,7 +79,7 @@ p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,
 t2 = seq(from=0,to=5000,by=1) 
 
 
-N0 = c(1,1,1000,1000,333,333,333,333,333,333)
+N0 = c(500,500,500,500,333,333,333,333,333,333)
 out2 = ode(y=N0,times=t2,func=SEAI,parms=p2)
 
 #plot(out2[,1],(out2[,2]+out2[,3]+out2[,4]+out2[,5]))
@@ -92,15 +92,15 @@ plot(out2[,1],(out2[,4]),type="l", xlab = "Time", ylab="Population", main = "b=0
 lines(out2[,1],(out2[,5]),col="RED")
 legend(x=40,y=6,legend=c("PRU","PRD"),bty="n",lwd = c(2,1),col=c("BLACK","RED"))
 
-plot(out2[,1],(out2[,6]),ylim=c(0,max(out2[,7])),type="l", xlab="Time", ylab="Population", main = "b=0.2239*365,d=0.30")
+plot(out2[,1],(out2[,6]),ylim=c(0,max(out2[,6])),type="l", xlab="Time", ylab="Population", main = "b=0.2239*365,d=0.30")
 lines(out2[,1],(out2[,7]),col="RED")
 legend(x=40,y=1000,legend=c("NU","ND"),bty="n",lwd=c(2,1),col=c("BLACK","RED"))
 
-plot(out2[,1],(out2[,8]),ylim=c(0,max(out2[,9])),type="l", xlab="Time",ylab="Population", main = "b=0.2239*365,d=0.30")
+plot(out2[,1],(out2[,8]),ylim=c(0,max(out2[,8])),type="l", xlab="Time",ylab="Population", main = "b=0.2239*365,d=0.30")
 lines(out2[,1],(out2[,9]),col="RED")
 legend(x=40,y=1000,legend=c("SU","SD"),bty="n",lwd=c(2,1),col=c("BLACK","RED"))
 
-plot(out2[,1],out2[,10],ylim=c(0,max(out2[,10])),type="l",xlab="Time",ylab="Asymptomatic Population (A)", main="b=0.2239*365,d=0.30")
+plot(out2[,1],out2[,10],ylim=c(0,max(out2[,11])),type="l",xlab="Time",ylab="Asymptomatic Population (A)", main="b=0.2239*365,d=0.30")
 lines(out2[,1],out2[,11],col="RED")
 
 plot(out2[,1],(out2[,8]+out2[,9])/(out2[,10]+out2[,11]))
@@ -116,13 +116,15 @@ bvalues = c(0,0.0001,0.0005,0.05,0.10,0.15,0.20,0.25,0.27,0.30,0.32,0.35,0.37,0.
 #bvalues = c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)
 wrvalues = c(0,0.10,0.15,0.20,0.25,0.27,0.30,0.32,0.35,0.37,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,1)
 dvalues = c(0,0.10,0.15,0.20,0.25,0.27,0.30,0.32,0.35,0.37,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,1)
+muSUvalues = c(0,0.10,0.15,0.20,0.25,0.27,0.30,0.32,0.35,0.37,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,1)
 for (i in 1:22){
-  b=bvalues[i]*1
+  b=bvalues[i]*0.01
   #wr=wrvalues[i]
-  #d=dvalues[i]
+  #d1=dvalues[i]
+  #muSU = muSUvalues[i]*0.01
   p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att)
   out3 = ode(y=N0,times=t2,func=SEAI,parms=p2)
-  out3eq[i] = out3[3999,2] + out3[3999,3]+out3[3999,4] + out3[3999,5]
+  out3eq[i] = out3[3999,3]+out3[3999,2] +out3[3999,4] + out3[3999,5]
   out3eq2[i] = (out3[3999,4]+out3[3999,5])/(out3[3999,2] + out3[3999,3] + out3[3999,4] + out3[3999,5])
   out3eq3[i] = out3[3999,6]+out3[3999,7]
   out3eq4[i] = out3[3999,8]+out3[3999,9]
@@ -131,9 +133,9 @@ for (i in 1:22){
 }
 
 #plot(bvalues,out3eq2,xlab="b/10",ylab="Proportion resistant")
-plot(bvalues,out3eq2)#,ylab="Proportion Resistant")
+plot(bvalues,out3eq)#,ylab="Proportion Resistant")
 #print(out3eq2)
-print(out3eq2)
+print(out3eq)
 
 mat6 = matrix(data=0,nrow=22,ncol=22)
 #mat7 = matrix(data=0,nrow=22,ncol=22)
@@ -145,8 +147,9 @@ for (i in 1:22){
   for (j in 1:22){
     #d = dvalues[i]
     #    A = Avalues[j]
-    wr = wrvalues[i]
-    b = bvalues[j]*1
+    #wr = wrvalues[i]
+    muSU = muSUvalues[i]*0.01
+    b = bvalues[j]*0.01
     p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d=d,m=m)
     out5 = ode(y=N0,times=t2,func=SEAI,parms=p2)
     #    PTotal = out5[3999,2] + out5[3999,3] + out5[3999,4] + out5[3999,5]