diff --git a/ODESSeasonality (1).R b/ODESSeasonality (1).R
new file mode 100644
index 0000000..4997f9d
--- /dev/null
+++ b/ODESSeasonality (1).R	
@@ -0,0 +1,454 @@
+require(deSolve)
+#library(tidyverse)
+SEAI<-function(t,y,p){
+  PSU = y[1]
+  PSD = y[2]
+  #PS = y[3]
+  PRU = y[3]
+  PRD = y[4]
+  #PR = y[6]
+  NU = y[5]
+  ND = y[6]
+  SU = y[7]
+  SD = y[8]
+  AU = y[9]
+  AD = y[10]
+  x = y[11]
+  beta=y[12]
+  with(as.list(p),{
+    
+    #gammainit = vector()
+    #gamma2=vector()
+    #gammainit1= 0.3438303
+    #numsteps = 10
+    #for (i in 1:numsteps){
+    #ratr1 = (PRU+PRD)/(PSU+PSD+PRU+PRD)
+    #P1 = ratr1
+    #TauRN1=((P1*((1-d1)+(d1*TauRD1)))-(d1*TauRD1))/(1-d1) #calculate TauRN1 from equation A4
+    #gammainit = gammainit1
+    #for(k in 1:cmax){
+    #  for(l in 1:hmax+1){
+    #    if(h1[l,k]>c1[k]){
+    #      x1[l,k]=0
+    #      Tau1[l,k]=0
+    #    }else if(h1[l,k]==c1[k]){
+    #      x1[l,k]=1
+    #      Tau1[l,k]=1
+    #    }else{
+    #      x1[l,k]=(h1[l,k]*(1-S))/(((c1[k]-h1[l,k]))+(h1[l,k]*(1-S)))
+    #      Tau1[l,k]>0
+    #      Tau1[l,k]=((x1[l,k])^2)+(2*x1[l,k]*(1-x1[l,k])*((1+gammainit)/2)^2)+(((1-x1[l,k])^2)*(1-(1#/(c1[k]-h1[l,k])))*(gammainit*gammainit))
+    #    }
+    #  }
+    #}
+    
+    
+    
+    
+    #n = 2
+    
+    #nprob = pbinom(h1[r,q],c1[q],prob=p)
+    #nprob = matrix(data=0,nrow=hmax+1,ncol=cmax)
+    #dn2 = ((1+gammainit)/2)^n
+    #Seq8=matrix(data=0,nrow=hmax+1,ncol=cmax)
+    #Seq9=matrix(data=0,nrow=hmax+1,ncol=cmax)
+    #for(q in 1:cmax){
+    #  for(r in 1:(hmax+1)){
+    #    if (h1[r,q]<=c1[q]){
+    #      Seq8[r,q]=(freqc1[q])
+    #    }else{
+    #      Seq8[r,q]=0
+    #    }
+    #    if(h1[r,q]<c1[q]){
+    #      nprob[r,q] = pbinom(h1[r,q],c1[q],prob=P1)
+    #      Seq9[r,q]=(freqc1[q]*nprob[r,q]*(1-(1/(c1[q]-h1[r,q])))*((1-x1[r,q])^2)*(gammainit^n)) #numerator of equation 1 with Tau=1 for all cases except h=0
+    #    }else{
+    #      Seq9[r,q]=0
+    #    }
+    #  }
+    #}
+    #alpha=sum(Seq9)/(sum(Seq8)) #Equation 1
+    #print(Seq7)
+    #print(alpha)
+    
+    #print("alpha is")
+    #print(alpha)
+    
+    #nprob2 = matrix(data=0,nrow=hmax+1,ncol=cmax)
+    #Seq10=matrix(data=0,nrow=hmax+1,ncol=cmax)
+    #Seq11=matrix(data=0,nrow=hmax+1,ncol=cmax)
+    #for(q in 1:cmax){
+    #  for(r in 1:(hmax+1)){    
+    #if (h1[r,q]<=c1[q]){
+    # Seq10[r,q]=(freqc1[q])
+    #}else{
+    #  Seq10[r,q]=0
+    #}
+    #if (h1[r,q]<=c1[q]){
+    #   nprob2[r,q] = pbinom(h1[r,q],c1[q],prob=P1)
+    #    Seq11[r,q]=(freqc1[q]*nprob2[r,q]*(1-dn2)*((1-x1[r,q])*x1[r,q])) #numerator of equation 1 with Tau=1 for all cases except h=0
+    #   }else{
+    #      Seq11[r,q]=0
+    #    }
+    #  }
+    #}
+    #beta= sum(Seq11)/(sum(Seq10))
+    
+    #print("beta is")
+    #print(beta)
+    
+    
+    #nprob3 = matrix(data=0,nrow=hmax+1,ncol=cmax)
+    #Seq12=matrix(data=0,nrow=hmax+1,ncol=cmax)
+    #Seq13=matrix(data=0,nrow=hmax+1,ncol=cmax)
+    #for(q in 1:cmax){
+    #  for(r in 1:(hmax+1)){ #denominator of Equation 1
+    #    if (h1[r,q]<=c1[q]){
+    #      Seq12[r,q]=(freqc1[q])
+    #    }else{
+    #      Seq12[r,q]=0
+    #    }
+    #    if (h1[r,q]<=c1[q]){
+    #      nprob3[r,q] = pbinom(h1[r,q],c1[q],prob=gammainit)
+    #      Seq13[r,q]=((freqc1[q]*nprob3[r,q]*h1[r,q]*(1-S))/(h1[r,q]*(1-S)+(c1[q]-h1[r,q]))) #numerator of equation 1 with Tau=1 for all cases except h=0
+    #    }else{
+    #      Seq13[r,q]=0
+    #    }
+    
+    #  }
+    #}
+    #Taubar= sum(Seq13)/sum(Seq12)#Equation 1
+    
+    #print("Taubar is")
+    #print(Taubar)
+    
+    #Seq14 = vector()
+    #Seq15 = matrix(data=0,nrow=hmax+1,ncol=cmax)
+    #for(q in 1:cmax){
+    # Seq14[q]=(freqc1[q]) #denominator of Equation 1
+    #  for(r in 1:(hmax+1)){
+    #    Seq15[r,q]=freqc1[q]*freqh1[r,q]*c1[q]#*Tau1[r,q]
+    #    
+    #  }
+    #}
+    
+    #TB = sum(Seq15)/sum(Seq14)
+    
+    
+    
+    #mu=1e-8
+    #gamma100 = ((gammainit*Taubar*(1-mu))+((1-gammainit)*(1-Taubar)*mu))/((gammainit*Taubar)+((1#-gammainit)*(1-Taubar)))
+    #print("gamma100 is")
+    #print(gamma100)
+    #pprime = (gamma100+beta-alpha)/(1+beta-alpha)
+    #print("pprime is")
+    #print(pprime)
+    #gamma= ratr #determined from popgen model
+    #wr[j] = 1-S[j]
+    #gammainit[i+1] = gamma2[i]
+    #}
+    #L = 10
+    #kap = 1
+    a =365/(2*pi)
+    ps =0
+    b = abs(beta)#+(bavgnull*L)+bavgnull)
+    #b =cos(kap*sin(t/a)-(t/a))-L
+    
+    
+    N = y[5]+y[6]+y[7]+y[8]+y[9]+y[10]
+    m = (PSU+PSD+PRU+PRD)/N
+    I = 1-exp(-m)#(1-((PSU+PSD+PRU+PRD)/(N*K))) 1-exp(-m)
+    ratr = ((PRU+PRD)*wr)/(PSU+PSD+((PRU+PRD)*wr))
+    ratr1 = (PRU+PRD)/(PSU+PSD+PRU+PRD)
+    rats = (PSU+PSD)/(PSU+PSD+((PRU+PRD)*wr))
+    rats1 = (PSU+PSD)/(PSU+PSD+PRU+PRD)
+    #P1 = ratr1
+    #ps = 0.34
+    #gammainit = (ratr1*1)+(rats*ps)
+    #gamma = (ratr1*1)+(rats*((gammainit - ratr1)/rats))
+    #A = ((1+gamma)/2)^2 #prob that recombinant is resistant given RxS recombination
+    #lambda = (PSU+PSD)/N
+    #cr = 1-((lambda*exp(-lambda))/(1-exp(-m)))#-exp(-lambda))#/(1-exp(-m)))
+    #B = (1-(1/cr))*(gamma^2) #prob that recombinant is resistant given SxS recombination
+    #rats1 = (PSU+PSD-1)/(PSU+PSD+((PRU+PRD)*wr))
+    #B = cr*(gamma^2)
+    #sigma1 = (ratr*A)+(rats*B)
+    #sigma2 = (rats*(1-A))
+    
+    eps = (PSD+PSU+((PRD+PRU)*(wr)))/(PSD+PSU+PRD+PRU)
+    r = 0.5*(1-(((m*exp(-m)))/(1-exp(-m))))
+    #r=0
+    upsilon = (1/K)*(K-((PSU+PSD+PRU+PRD)/N))
+    delR = b*upsilon*(PRU+PRD)*wr#*I
+    delS = b*upsilon*(PSU+PSD)*ws#*I
+    delP = delS +delR
+    n=2
+    #TauRN1=((P1*((1-d1)+(d1*TauRD1)))-(d1*TauRD1))/(1-d1) #calculate TauRN1 from equation A4
+    #FUNC2<-function(smallp,E,S,Tau1,x1,f1,c1,h1,freqc1,freqh1){
+    # for(k in 1:cmax){
+    #  for(l in 1:hmax+1){
+    #   if(h1[l,k]>c1[k]){
+    #    x1[l,k]=0
+    #   Tau1[l,k]=0
+    #}else if(h1[l,k]==c1[k]){
+    # x1[l,k]=1
+    #Tau1[l,k]=1
+    #}else{
+    # x1[l,k]=(h1[l,k]*(1-S))/(((c1[k]-h1[l,k]))+(h1[l,k]*(1-S)))
+    #Tau1[l,k]>0
+    #Tau1[l,k]=((x1[l,k])^2)+(2*x1[l,k]*(1-x1[l,k])*((1+smallp)/2)^2)+(((1-x1[l,k])^2)*(1-(1/(c1[k]-h1[l,k])))*(smallp*smallp))
+    #}
+    #}
+    #}
+    #return(((freqc1[1]*((freqh1[1,1]*Tau1[1,1])+(freqh1[2,1]*Tau1[2,1])))+(freqc1[2]*((freqh1[1,2]*Tau1[1,2])+(freqh1[2,2]*Tau1[2,2])+(freqh1[3,2]*Tau1[3,2])))+(freqc1[3]*((freqh1[1,3]*Tau1[1,3])+(freqh1[2,3]*Tau1[2,3])+(freqh1[3,3]*Tau1[3,3])+(freqh1[4,3]*Tau1[4,3])))+(freqc1[4]*((freqh1[1,4]*Tau1[1,4])+(freqh1[2,4]*Tau1[2,4])+(freqh1[3,4]*Tau1[3,4])+(freqh1[4,4]*Tau1[4,4])+(freqh1[5,4]*Tau1[5,4])))+(freqc1[5]*((freqh1[1,5]*Tau1[1,5])+(freqh1[2,5]*Tau1[2,5])+(freqh1[3,5]*Tau1[3,5])+(freqh1[4,5]*Tau1[4,5])+(freqh1[5,5]*Tau1[5,5])+(freqh1[6,5]*Tau1[6,5])))+(freqc1[6]*((freqh1[1,6]*Tau1[1,6])+(freqh1[2,6]*Tau1[2,6])+(freqh1[3,6]*Tau1[3,6])+(freqh1[4,6]*Tau1[4,6])+(freqh1[5,6]*Tau1[5,6])+(freqh1[6,6]*Tau1[6,6])+(freqh1[7,6]*Tau1[7,6]))+(freqc1[7]*((freqh1[1,7]*Tau1[1,7])+(freqh1[2,7]*Tau1[2,7])+(freqh1[3,7]*Tau1[3,7])+(freqh1[4,7]*Tau1[4,7])+(freqh1[5,7]*Tau1[5,7])+(freqh1[6,7]*Tau1[6,7])+(freqh1[7,7]*Tau1[7,7])+(freqh1[8,7]*Tau1[8,7]))+(freqc1[8]*((freqh1[1,8]*Tau1[1,8])+(freqh1[2,8]*Tau1[2,8])+(freqh1[3,8]*Tau1[3,8])+(freqh1[4,8]*Tau1[4,8])+(freqh1[5,8]*Tau1[5,8])+(freqh1[6,8]*Tau1[6,8])+(freqh1[7,8]*Tau1[7,8])+(freqh1[8,8]*Tau1[8,8])+(freqh1[9,8]*Tau1[9,8])))))/((freqc1[1])+(freqc1[2])+(freqc1[3])+(freqc1[4])+(freqc1[5])+(freqc1[6])+(freqc1[7])+(freqc1[8])))-TauRN1)
+    #}
+    #Recomb = uniroot(f=FUNC2,interval=c(0,1),S=S,Tau1=Tau1,x1=x1,c1=c1,h1=h1,freqc1=freqc1,freqh1=freqh1,f1=f1) #find root of equation 1 between 0 and 1
+    #gammainit=(Recomb$root)
+    
+    #gammainit3 = 0.471481
+    #gamma2final = gamma2[numsteps]
+    #PSU to PRU: (sigma1*b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
+    #PSU to PRD: (sigma1*b*(PSU/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
+    #PSD to PRU: (sigma1*b*(PSD/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
+    #PSD to PRD: (sigma1*b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
+    #PRU to PSU: (sigma2*b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
+    #PRU to PSD: (sigma2*b*(PRU/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
+    #PRD to PSU: (sigma2*b*(PRD/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
+    #PRD to PSD: (sigma2*b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
+    #PSD to PSD: ((1-sigma1)*b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
+    #PSU to PSU: ((1-sigma1)*b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
+    #PRD to PRD: ((1-sigma2)*b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
+    #PRU to PRU: ((1-sigma2)*b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))
+    #dPSU.dt = -((1-sigma1)*b*((PSU+PSD)/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRD/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (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)#-(del*att*(1/N)*(NU+SU+AU)*PSU)
+    #dPSD.dt = -((1-sigma1)*b*((PSU+PSD)/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma2*b*(PRU/N)*wr*(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)) - (b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSD)*(PSD)#-(del*att*(1/N)*(ND+SD+AD)*PSD)
+    #dPRU.dt = -((1-sigma2)*b*((PRU+PRD)/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSD/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))- (b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRU)*(PRU)#-(del*att*(1/N)*(NU+SU+AU)*PRU)
+    #dPRD.dt = -((1-sigma2)*b*((PRU+PRD)/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSU/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))- (b*(PRU/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRD)*(PRD)#-(del*att*(1/N)*(ND+SD+AD)*PRD)
+    
+    gammavec = c((n/n-1),x)
+    gammavec2 = c(0,x)
+    if (x<0){
+      gamma = max(gammavec2)
+    }else{
+      gamma = min(gammavec)
+    }
+    #x1 = ((PRU+PRD+PSU+PSD)*x-(PRU+PRD)*(1-wr))/(PSD+PSU+((PRD+PRU)*(wr)))
+    #x1 = (x2^n)/(1+x2^n)
+    x1 = gamma*(rats1)+ratr1*wr
+    dPSU.dt = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSU*PSU
+    dPSD.dt = ((ND+SD+AD)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSD*PSD
+    dPRU.dt = ((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x1^n))-muRU*PRU
+    dPRD.dt = ((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x1^n))-muRD*PRD
+    e = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSU*PSU+((ND+SD+AD)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSD*PSD+((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x1^n))-muRU*PRU+((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x1^n))-muRD*PRD
+    er = ((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x1^n))-muRU*PRU+((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x1^n))-muRD*PRD
+    es = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSU*PSU+((ND+SD+AD)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSD*PSD
+    #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*ratr*ND*I*b*eps) - (del*att*ND)
+    dSU.dt = (b*I*eps*NU*rho1)+(SD/Tau)-(b*I*eps*SU*(d2+rho3)) - (del*att*SU)
+    dSD.dt = (ratr*b*I*eps*ND*rho2)+(b*I*eps*SU*d2)-(SD/Tau)-(ratr*b*I*eps*SD*rho4) - (del*att*SD)
+    #dAU.dt = rat*(AU+AD)*(1-((AU+AD)/L))+(AD/Tau)+(b*I*eps*SU*rho3)-(b*I*eps*omega1*AU*d) - (del*att*AU)
+    dAU.dt = (AD/Tau)+(b*I*eps*SU*rho3)-(b*I*eps*omega1*AU*d3) - (del*att*AU)
+    #dAD.dt = rat*(AU+AD)*(1-((AU+AD)/L))+(b*I*eps*SD*rho4)+(b*I*eps*omega1*AU*d)-(AD/Tau) - (del*att*AD)
+    dAD.dt = (ratr*b*I*eps*SD*rho4)+(b*I*eps*omega1*AU*d3)-(AD/Tau) - (del*att*AD)
+    #dgamma.dt = (1/(PSU+PSD+((1-S)*(PRU+PRD))))*((1/rats)-(gamma))*(dPRD.dt+dPRU.dt)
+    #dgamma.dt = ratr + rats*((gamma)-ratr)-((1/(PSU+PSD+((1-S)*(PRU+PRD))))*(dPRD.dt+dPRU.dt)*((gamma)-ratr))
+    #dgamma.dt = -(dPRD.dt+dPRU.dt)/(dPSD.dt+dPSU.dt)
+    #dx.dt = (1/((PSU+PSD+PRU+PRD)+(e)))*((er)-(x*(e))+((es)*((x*(PSD+PSD+PRU+PRD)-(PRU+PRD))/(PSU+PSD))))
+    dx.dt = (delP*(x1-x)+delR*(x-1))/((PSU+PSD)+delS)
+    dbeta.dt = ((L*bavgnull*sin(((t/a)+ps)-kap*sin((t/a)+ps))*(kap*cos((t/a)+ps)-1))/a)
+    return(list(c(dPSU.dt,dPSD.dt,dPRU.dt,dPRD.dt,dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,dx.dt,dbeta.dt)))
+  })
+}
+
+#b = 10
+
+nstrains = 150
+
+del = 1
+att = 1/365/30#.000000000001
+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
+rho3 = 1/nstrains #prob of immunity gain from SU to AU
+rho4 = rho3#=1/nstrains #prob of immunity gain from SD to AD
+omega1 = 5/100 #prob of appearance of new antigen
+omega2 = omega1 #prob of appearance of new antigen only in resistant strains
+
+bavgnull = 0.1
+L =  0.8
+kap = -1
+
+ws = 1
+wr = 1-0.03
+S = 1-wr
+#muSU = 1
+muSU=1/250#0.001
+#muSD = 20
+muSD=0.99
+#muRD = 1
+muRD = muSU#0.001
+#muRU = 1
+muRU = muSU#0.001
+Tau = 20#/365
+d1 = 0.87#*0.3
+d2 = d1#/6
+d3 = d1#/10
+#m = 3
+K=10
+p2 = list(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, K=K, bavgnull=bavgnull)
+
+tmax=40000
+t2 = seq(from=0,to=tmax,by=1)
+
+
+N0 = c(10000,10000,1,1,333,333,333,333,333,333,0.34,bavgnull+(bavgnull*L))
+out2 = ode(y=N0,times=t2,func=SEAI,parms=p2)
+
+plot(out2[,1],(out2[,2]),ylim=c(0,max(out2[,2])),type="l", xlab="Time", ylab="Population", main = paste("b=",b,", d=",d1))
+lines(out2[,1],(out2[,3]),col="RED")
+legend("bottomright",legend=c("PSU","PSD"),bty="n",lwd = c(2,1),col=c("BLACK","RED"))
+
+plot(out2[,1],(out2[,4]),type="l", ylim=c(0,max(out2[,5])),xlab = "Time", ylab="Population", main = "b=0.2239*365,d=0.30")
+lines(out2[,1],(out2[,5]),col="RED")
+legend("bottomright",legend=c("PRU","PRD"),bty="n",lwd = c(2,1),col=c("BLACK","RED"))
+
+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("bottomright",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")
+lines(out2[,1],(out2[,9]),col="RED")
+legend("bottomright",legend=c("SU","SD"),bty="n",lwd=c(2,1),col=c("BLACK","RED"))
+
+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[,12],type="l")
+#print(out2[,12][35000])
+
+plot(out2[,1],abs(out2[,13]),type="l",xlim=c(0,2000))
+
+
+plot(out2[,1],(out2[,4]+out2[,5])/(out2[,2]+out2[,3]+out2[,4]+out2[,5]),type="l",xlim=c(30000,40000),ylim=c(0.75,0.85))
+
+#out3eq = vector(mode = "numeric", length = 18)
+out3eq2 = vector(mode = "numeric", length = 18)
+#out3eq3 = vector(mode = "numeric", length = 18)
+#out3eq4 = vector(mode = "numeric", length = 18)
+#out3eq5 = vector(mode="numeric", length=18)
+#out3eq7 = vector(mode="numeric", length=18)
+#out4 = 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(-3,-2.5,-2,-1.5, -1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2,0.4,0.6,0.8,1,1.5,2)#,3,4)
+bavg = vector(mode="numeric",length=13)
+#bvalues = c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2)
+kapvalues = c(-1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 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)
+#bvalues = seq(from=-2.5,to=0.5,by=1)
+wrvalues = c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)
+#wrvalues = c(0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1)
+#wrvalues = c(0.20,0.22,0.24,0.26,0.28,0.30,0.32,0.34,0.36,0.38,0.40)
+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)
+delvalues = c(10,20,40,50,60,70,90,100,150,200,250,300,350,400,450,500,550,600,700,800,900,1000)
+#rho3values = c(0.0005,0.002,0.004,0.005,0.007,0.008,0.01,0.02,0.04,0.05,0.06,0.07,0.08,0.09,0.10,0.11)
+rho3values = c(0.0005,0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01)
+omegavalues = c(0.01,0.012,0.014,0.016,0.018,0.02,0.022,0.024,0.026,0.028,0.030,0.032,0.034,0.036,0.038,0.040)#,0.042,0.044,0.046,0.048,0.050,0.052)
+#Lvalues = c(-3,-2.5,-2,-1.5, -1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2,0.4,0.6,0.8,1,1.5,2)#,3,4)
+Lvalues = c(0.04,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)
+bavgvalues = c(-1.5, -1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.5,2)
+for (i in 1:18){
+  #del=delvales[i]
+  #b=10^(bvalues[i])
+  #wr=wrvalues[i]
+  #d1=dvalues[i]
+  #d2=d1
+  #d3=d1
+  #muSU = muSUvalues[i]*0.1
+  #rho3 = rho3values[i]
+  #rho4=rho3
+  bavgnull = 10^bavgvalues[i]
+  #bavg[i] = bavgnull*L
+  #L=Lvalues[i]
+  #kap=kapvalues[i]
+  #omega1 = omegavalues[i]
+  p2 = list(b=b,ws=ws,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,kap=kap,bavgnull=bavgnull)
+  out3 = ode(y=N0,times=t2,func=SEAI,parms=p2)
+  #out3eq[i] = out3[39999,3]+out3[39999,2] +out3[39999,4] + out3[39999,5]
+  out3eq2[i] = (out3[39999,4]+out3[39999,5])/(out3[39999,2] + out3[39999,3] + out3[39999,4] + out3[39999,5])
+  #out3eq3[i] = (out3[3999,6]+out3[3999,7])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11])
+  #out3eq4[i] = (out3[39999,8]+out3[39999,9])/(out3[39999,6]+out3[39999,7]+out3[39999,8]+out3[39999,9]+out3[39999,10]+out3[39999,11])
+  #out3eq5[i] =(out3[39999,10]+out3[39999,11])/(out3[39999,6]+out3[39999,7]+out3[39999,8]+out3[39999,9]+out3[39999,10]+out3[39999,11])
+  #out4 [i] = (1/(K*(out3[399999,6]+out3[399999,7]+out3[399999,8]+out3[399999,9]+out3[399999,10]+out3[399999,11])))*(-(((out3[3999,3]+out3[3999,2]+out3[3999,4] + out3[3999,5])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11]))-K))
+  #out4[i] =  (1/(K*(K)))*(-(((out3[3999,3]+out3[3999,2]+out3[3999,4] + out3[3999,5])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11]))-K))
+  #out3eq6[i] = (out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11])
+  #out3eq7[i] = (out3[39999,7]+out3[39999,9]+out3[39999,11])/(out3[39999,6]+out3[39999,7]+out3[39999,8]+out3[39999,9]+out3[39999,10]+out3[39999,11])
+}
+
+#plot(bvalues,out3eq2,xlab="b/10",ylab="Proportion resistant")
+#plot(bvalues,1-out4)#,ylab="Proportion Resistant")
+plot(bavgvalues,out3eq2)
+#print(out3eq2)
+print(out3eq2)
+print(bavg)
+
+
+mat6 = matrix(data=0,nrow=11,ncol=18)
+#mat7 = matrix(data=0,nrow=11,ncol=17)
+#mat8 = matrix(data=0,nrow=11,ncol=19)
+#mat9 = matrix(data=0,nrow=11,ncol=18)
+#mat10 = matrix(data=0,nrow=11,ncol=18)
+#mat11 = matrix(data=0,nrow=11,ncol=18)
+#mat12 = matrix(data=0,nrow=11,ncol=18)
+for (i in 1:11){
+  for (j in 1:18){
+    #d = dvalues[i]
+    #    A = Avalues[j]
+    #rho3=rho3values[i]
+    #rho4=rho3
+    #omega1 = omegavalues[i]
+    bavgnull = 10^bavgvalues[j]
+    wr = wrvalues[i]
+    #kap=kapvalues[i]#muSU = muSUvalues[i]*0.01
+    #b = 10^(bvalues[j])
+    #L=Lvalues[j]
+    #omega1 = omegavalues[j]
+    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,bavgnull=bavgnull)
+    out5 = ode(y=N0,times=t2,func=SEAI,parms=p2)
+    #    PTotal = out5[3999,2] + out5[3999,3] + out5[3999,4] + out5[3999,5]
+    #    mat[i,j] = (out5[3999,4] + out5[3999,5])/PTotal
+    #mat2[i,j] = (out5[3999,2] + out5[3999,3] + out5[3999,4] + out5[3999,5])/6000
+    #mat3[i,j] = out5[3999,6]/out5[3999,7]
+    mat6[i,j] = (out5[39999,4]+out5[39999,5])/(out5[39999,2] + out5[39999,3] + out5[39999,4] + out5[39999,5])
+    #mat7[i,j] = out5[3999,2]+out5[3999,3]+out5[3999,4]+out5[3999,5]
+    #mat8[i,j]=(out5[39999,10]+out5[39999,11])/(out5[39999,6]+out5[39999,7]+out5[39999,8]+out5[39999,9]+out5[39999,10]+out5[39999,11])
+    #mat9[i,j] = (out5[3999,8]+out5[3999,9])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11])
+    #mat10[i,j] = (out5[3999,6]+out5[3999,7])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11])
+    #mat11[i,j] = (out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11])
+    #mat12[i,j]=(out5[3999,7]+out5[3999,9]+out5[3999,11])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11])
+  }
+}
+print(mat6)
+#print(mat7)
+#print(mat8)
+#print(mat9)
+#print(mat10)
+#print(mat11)
+#print(mat12)
+
+
+FLUC<-function(t,L,kap,beta,a){
+  L = 1
+  beta = 10^-0.5
+  kap = 1
+  a =1/(2*pi)
+  b = -L*beta*cos(kap*sin(t/a)-(t/a))+beta*2
+  return(b)
+}
+
+integrate(FLUC,0,1)
+
+
+
+