diff --git a/ODEsImmunity.R b/ODEsImmunity.R
index 4ec90dc..27568d5 100644
--- a/ODEsImmunity.R
+++ b/ODEsImmunity.R
@@ -15,28 +15,31 @@ SEAI<-function(t,y,p){
   AD = y[10]
   with(as.list(p),{
     N = y[5]+y[6]+y[7]+y[8]+y[9]+y[10]
+    m = (PSU+PSD+PRU+PRD)/N
+    I = 1-exp(-m)
     eps = (((PSD+PSU+((PRD+PRU)*(wr))))/(PSD+PSU+PRD+PRU))
     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)
-    dNU.dt = (del*N*(1-(N/L))) + (ND/Tau)-(b*eps*NU*(d+rho1))-(muNU*NU)# - (del*NU)
-    dND.dt = (b*eps*NU*d)-(ND/Tau)-(b*eps*ND*rho2)-(muND*ND)# - (del*ND)
-    dSU.dt = (b*eps*omega1*AU*SU)+(b*eps*NU*rho1)+(SD/Tau)-(b*eps*SU*(d+rho3))# - (del*SU)
-    dSD.dt = (b*eps*omega2*AD*SD)+(b*eps*ND*rho2)+(b*eps*SU*d)-(SD/Tau)-(b*eps*SD*rho4)# - (del*SD)
-    dAU.dt = (AD/Tau)+(b*eps*SU*rho3)-(b*eps*omega1*AU*SU)# - (del*AU)
-    dAD.dt = (b*eps*SD*rho4)-(b*eps*omega2*AD*SD)-(AD/Tau)# - (del*AD)
+    dNU.dt = (del*N)+ (ND/Tau)-(b*I*eps*NU*(d+rho1))-(muNU*NU) - (del*att*NU)
+    dND.dt = (b*I*eps*NU*d)-(ND/Tau)-(b*I*eps*ND*rho2)-(muND*ND) - (del*att*ND)
+    dSU.dt = (b*I*eps*NU*rho1)+(SD/Tau)-(b*I*eps*SU*(d+rho3)) - (del*att*SU)
+    dSD.dt = (b*I*eps*ND*rho2)+(b*I*eps*SU*d)-(SD/Tau)-(b*I*eps*SD*rho4) - (del*att*SD)
+    dAU.dt = (AD/Tau)+(b*I*eps*SU*rho3)-(b*I*eps*omega1*AU*d) - (del*att*AU)
+    dAD.dt = (b*I*eps*SD*rho4)+(b*I*eps*omega1*AU*d)-(AD/Tau) - (del*att*AD)
     return(list(c(dPSU.dt,dPSD.dt,dPRU.dt,dPRD.dt,dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt)))
   })
 }
 
-b = 0.2239*365
+b = 20000
 
 nstrains = 20
 
-del = (1/60)
-muND = 0.05*20
-muNU = 0.1*20
+att = 0
+del = 0
+muND = 0
+muNU = 0
 rho1 = 0.5 #prob of immunity gain from NU to SU
 rho2 = 0.5 #prob of immunity gain from ND to SD
 rho3 = 1/nstrains #prob of immunity gain from SU to AU
@@ -45,32 +48,32 @@ omega1 = 1/100 #prob of appearance of new antigen
 omega2 = 1/100 #prob of appearance of new antigen only in resistant strains
 
 
-L = 5000
+L = 10000
 
 ws = 1
 wr = 1-0.6
-muSU = 1
-#muSU=0.05
-muSD = 20
-#muSD=0.99
-muRD = 1
-#muRD = 0.05
-muRU = 1
-#muRU = 0.05
-Tau = 20/365
+#muSU = 1
+muSU=0.05
+#muSD = 20
+muSD=0.99
+#muRD = 1
+muRD = 0.05
+#muRU = 1
+muRU = 0.05
+Tau = 20#/365
 d = 0.9
 #m = 3
-K=6
-p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d=d,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L)
+K=10
+p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d=d,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att)
 
 
-t2 = seq(from=0,to=1000,by=1/365) 
+t2 = seq(from=0,to=5000,by=1) 
 
 
-N0 = c(500,500,500,500,1000,1000,1000,1000,1000,1000)
+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])/6000)
+#plot(out2[,1],(out2[,2]+out2[,3]+out2[,4]+out2[,5]))
 
 plot(out2[,1],(out2[,2]),type="l", xlab="Time", ylab="Population", main = "b=0.2239*365,d=0.30")
 lines(out2[,1],(out2[,3]),col="RED")
@@ -80,35 +83,35 @@ 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]),type="l",ylim=c(0,max(out2[,6])), xlab="Time", ylab="Population", main = "b=0.2239*365,d=0.30")
+plot(out2[,1],(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]),type="l",ylim=c(0,max(out2[,8])), xlab="Time", ylab="Population", main = "b=0.2239*365,d=0.30")
+plot(out2[,1],(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],type="l",ylim=c(0,max(out2[,10])),xlab="Time",ylab="Asymptomatic Population (A)", main="b=0.2239*365,d=0.30")
+plot(out2[,1],out2[,10],type="l",xlab="Time",ylab="Asymptomatic Population (A)", main="b=0.2239*365,d=0.30")
 lines(out2[,1],out2[,11],col="RED")
 
 
-out3eq = vector(mode = "numeric", length = 21)
-out3eq2 = vector(mode = "numeric", length = 21)
-bvalues = c(0,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.70,0.80,0.90,1,2,3,4,5,6,7,8)
+out3eq = vector(mode = "numeric", length = 22)
+out3eq2 = vector(mode = "numeric", length = 22)
+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.40,0.45,0.50,0.55,0.60,0.70,0.80,0.90,1)
 #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)
 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)
-for (i in 1:21){
-  #b=bvalues[i]*365
+for (i in 1:22){
+  b=bvalues[i]*10
   #wr=wrvalues[i]
-  d=dvalues[i]
-  p3 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d=d,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L)
+  #d=dvalues[i]
+  p3 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d=d,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=p3)
   #out3eq[i] = out3[3999,2] + out3[3999,3]+out3[3999,4] + out3[3999,5]
-  out3eq2[i] = (out3[30000,4]+out3[30000,5])/(out3[30000,2] + out3[30000,3] + out3[30000,4] + out3[30000,5])
+  out3eq2[i] = (out3[3999,4]+out3[3999,5])/(out3[3999,2] + out3[3999,3] + out3[3999,4] + out3[3999,5])
 }
 
-plot(dvalues,out3eq2)
+plot(bvalues,out3eq2,xlab="b/10")
 print(out3eq2)