diff --git a/batchRuns/ODESImmunity_v2.R b/batchRuns/ODESImmunity_v2.R
new file mode 100644
index 0000000..9c1e452
--- /dev/null
+++ b/batchRuns/ODESImmunity_v2.R
@@ -0,0 +1,437 @@
+require(deSolve)
+library(tidyverse)
+args <- commandArgs(trailingOnly = TRUE)
+print(args)
+paramFile<-args[1]
+No<-as.integer(args[2])
+resultsFolder <-args[3]
+WildStrainOnly<-as.integer(args[4])
+#paramFile<-"paramList.csv"
+#No<-1
+#resultsFolder<-"./"
+
+SEAI<-function(t,y,p){
+
+  with(as.list(c(y,p)),{
+    # total host population size
+    N = NU+ND+SU+SD+AU+AD
+    # total parasite population size
+    P = PWU+PWD+PRU+PRD
+    
+    # parasite/host ratio
+    m = P/N
+    # prevalence
+    I = 1-exp(-m)
+    # converted cost of resistance considering competition
+    h = 1-(ResCost/(1+0.1*exp(-competition*(I-1))))
+    #adjusted parasite population, transmissibility lower in PR
+    PADJ = PWU+PWD+(PRU+PRD)*h
+    
+    # pop size +1 is to prevent N = 0 for undefined scenario
+      fracN = P_acc_N/(NU+ND+1)
+      fracS = P_acc_S/(SU+SD+1)
+      fracA = P_acc_A/(AU+AD+1)
+    
+    # probability that the infected strain has not been seen by the adaptive immune system
+    infN = exp(-fracN/nstrains)
+    infS = exp(-fracS/nstrains)
+    infA = exp(-fracA/nstrains)
+
+    #transmissibility is reduced to the number of infected hosts, 
+    #receiving hosts needs to have open "resource room" for accepting new strains
+    B = b*(I/m)*(1-m/capacity) 
+    
+    # per capita biting rate for hosts who are not treated
+    Bu = B*PADJ/N
+    # per capita biting rate for hosts who are drug treated
+    Bd = B*(PRU+PRD)*h/N
+    #print(B)
+    #eq b
+    infD = ND*infN+SD*infS+AD*infA
+    infU = NU*infN+SU*infS+AU*infA
+    #print(infD)
+    #print(infU)
+    
+    eq1 = B*h*(infU/N)
+    eq2 = B*h*(infD/N)
+    eq3 = B*(infU/N)
+    eq4 = B*(infD/N)
+    
+    # Parasite pop ODE
+    dPWU.dt = PWU*eq3 + PWD*eq3 - muWU*PWU
+    dPWD.dt = PWD*eq4 + PWU*eq4 - muWD*PWD
+    dPRU.dt = PRU*eq1 + PRD*eq1 - muRU*PRU
+    dPRD.dt = PRD*eq2 + PRU*eq2 - muRD*PRD
+    #sprintf("%.3f, %.3f,%.3f,%.3f",dPWU.dt,dPWD.dt,dPRU.dt,dPRD.dt)
+
+    # Generalized immunity ODE
+    dNU.dt = birth + 
+      (ND/Tau)- # drug loose effectiveness
+      (NU*Bu*infN*d1) -  #naive class getting infected, and treated
+      (NU*Bu*rho1) - #naive class getting enough # of infections to move to S
+      (NU*Bu*muNU)- #naive class die from infections
+      (NU*att)#normal death rate
+
+   
+    #print(dNU.dt)
+    
+    dND.dt = (NU*Bu*infN*d1) - #naive class getting infected, and treated
+    (ND/Tau) - # drug loose effectiveness
+    (ND*Bd*rho1) -#naive class getting enough # of infections to move to S
+    (ND*Bd*muND) -#naive class die from infections
+    (ND*att) #normal death rate
+    #print(dND.dt)
+    
+    dSU.dt =  (NU*Bu*rho1) +#naive class getting enough # of infections to move to S
+    (SD/Tau)-  # drug loose effectiveness
+    (SU*Bu*infS*d2)- #symptomatic class getting infected, show symptoms, and treated
+    (SU*Bu*rho2) -#symptomatic class getting enough # of infections to move to AS
+    (SU*att) #normal death rate
+    
+    dSD.dt = (ND*Bd*rho1)+ #naive class getting enough # of infections to move to S
+     (SU*Bu*infS*d2)- #symptomatic class getting infected, show symptoms, and treated
+     (SD/Tau) - # drug loose effectiveness
+     (SD*Bd*rho2)- #symptomatic class getting enough # of infections to move to AS
+     (SD*att) #normal death rate
+    
+    dAU.dt = (SU*Bu*rho2)+ #symptomatic class getting enough # of infections to move to AS
+      (AD/Tau)-  # drug loose effectiveness
+     (AU*Bu*infA*omega*d3)- #asymptomatic class getting infected, show symptoms, and treated
+     (AU*att) #normal death rate
+    
+    dAD.dt = (SD*Bd*rho2) + #symptomatic class getting enough # of infections to move to AS
+     (AU*Bu*infA*omega*d3)- #asymptomatic class getting infected, show symptoms, and treated
+     (AD/Tau)-  # drug loose effectiveness
+     (AD*att) #normal death rate
+    
+    #sprintf("%.0f, %.0f,%.0f,%.0f,%.0f,%.0f",dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt)
+    
+    # track the mean number of times infected per immunity class
+    dP_acc_N.dt = NU*Bu + ND*Bd - # total number of infection received
+     (NU*(Bu*muNU+att))*P_acc_N/(NU+ND+1)- #when naive U class die how many infections deducted
+     (ND*(Bd*muND+att))*P_acc_N/(NU+ND+1)- #when naive D class die many infections deducted
+     ((NU*Bu*rho1)+(ND*Bd*rho1))*P_acc_N/(NU+ND+1) #when naive class leave to S
+    
+    dP_acc_S.dt = ((NU*Bu*rho1)+(ND*Bd*rho1))*P_acc_N/(NU+ND+1)+ #when naive class leave to S
+     SU*Bu + SD*Bd- #total number of new infections
+     (SU+SD)*att*P_acc_S/(SU+SD+1)- # normal death
+     ((SD*Bd*rho2)+ (SU*Bu*rho2))*P_acc_S/(SU+SD+1) #Symptomatic move to Asymptomatic
+    
+    dP_acc_A.dt = ((SD*Bd*rho2)+ (SU*Bu*rho2))*P_acc_S/(SU+SD+1)+ #Symptomatic move to Asymptomatic
+     AU*Bu + AD*Bd- #total number of new infections
+     (AU+AD)*att*P_acc_A/(AU+AD+1) # normal death
+    
+    a <- 365/(2*pi)
+    db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a))
+    #db.dt =0
+    #sprintf("%.0f, %.0f,%.0f,%.0f",dP_acc_N.dt,dP_acc_S.dt,dP_acc_A.dt,dbeta.dt)
+
+    return(list(c(dPWU.dt,dPWD.dt,dPRU.dt,dPRD.dt,dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt)))
+  })
+}
+
+SEAI_p<-function(t,y,p){
+  
+  with(as.list(c(y,p)),{
+    # total host population size
+    N = NU+ND+SU+SD+AU+AD
+    # total parasite population size
+    #PW = PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD
+    #PR = PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD
+    #P = PW+PR
+    
+    # parasite/host ratio
+    rNU = (PW_NU+PR_NU)/(NU+1)
+    rND = (PW_ND+PR_ND)/(ND+1)
+    rSU = (PW_SU+PR_SU)/(SU+1)
+    rSD = (PW_SD+PR_SD)/(SD+1)
+    rAU = (PW_AU+PR_AU)/(AU+1)
+    rAD = (PW_AD+PW_AD)/(AD+1)
+    
+    # prevalence
+    I_NU = 1-exp(-rNU)
+    I_ND = 1-exp(-rND)
+    I_SU = 1-exp(-rSU)
+    I_SD = 1-exp(-rSD)
+    I_AU = 1-exp(-rAU)
+    I_AD = 1-exp(-rAD)
+    
+    # individual contribution from W or R
+    ic_NU <- c(exp(-rNU*PW_NU/(PW_NU+PR_NU+1)), exp(-rNU*PR_NU/(PW_NU+PR_NU+1)))
+    ic_ND <- c(exp(-rND*PW_ND/(PW_ND+PR_ND+1)), exp(-rND*PR_ND/(PW_ND+PR_ND+1)))
+    ic_SU <- c(exp(-rSU*PW_SU/(PW_SU+PR_SU+1)), exp(-rSU*PR_SU/(PW_SU+PR_SU+1)))
+    ic_SD <- c(exp(-rSD*PW_SD/(PW_SD+PR_SD+1)), exp(-rSD*PR_SD/(PW_SD+PR_SD+1)))
+    ic_AU <- c(exp(-rAU*PW_AU/(PW_AU+PR_AU+1)), exp(-rAU*PR_AU/(PW_AU+PR_AU+1)))
+    ic_AD <- c(exp(-rAD*PW_AD/(PW_AD+PR_AD+1)), exp(-rAD*PR_AD/(PW_AD+PR_AD+1)))
+    
+    # converted contribution of W or R with mixed infection cost
+    h_NU = c((1-ic_NU[1])*(exp(-rNU)/ic_NU[1]+(1-exp(-rNU)/ic_NU[1])*(1+ResCost/2)),(1-ic_NU[2])*(exp(-rNU)/ic_NU[2]*cloneCost+(1-exp(-rNU)/ic_NU[2])*(1-mixCost/2)))
+    h_ND = c((1-ic_ND[1])*(exp(-rND)/ic_ND[1]+(1-exp(-rND)/ic_ND[1])*(1+ResCost/2)),(1-ic_ND[2])*(exp(-rND)/ic_ND[2]*cloneCost+(1-exp(-rND)/ic_ND[2])*(1-mixCost/2)))
+    h_SU = c((1-ic_SU[1])*(exp(-rSU)/ic_SU[1]+(1-exp(-rSU)/ic_SU[1])*(1+ResCost/2)),(1-ic_SU[2])*(exp(-rSU)/ic_SU[2]*cloneCost+(1-exp(-rSU)/ic_SU[2])*(1-mixCost/2)))
+    h_SD = c((1-ic_SD[1])*(exp(-rSD)/ic_SD[1]+(1-exp(-rSD)/ic_SD[1])*(1+ResCost/2)),(1-ic_SD[2])*(exp(-rSD)/ic_SD[2]*cloneCost+(1-exp(-rSD)/ic_SD[2])*(1-mixCost/2)))
+    h_AU = c((1-ic_AU[1])*(exp(-rAU)/ic_AU[1]+(1-exp(-rAU)/ic_AU[1])*(1+ResCost/2)),(1-ic_AU[2])*(exp(-rAU)/ic_AU[2]*cloneCost+(1-exp(-rAU)/ic_AU[2])*(1-mixCost/2)))
+    h_AD = c((1-ic_AD[1])*(exp(-rAD)/ic_AD[1]+(1-exp(-rAD)/ic_AD[1])*(1+ResCost/2)),(1-ic_AD[2])*(exp(-rAD)/ic_AD[2]*cloneCost+(1-exp(-rAD)/ic_AD[2])*(1-mixCost/2)))
+    
+    #transmissibility is reduced to the number of infected hosts, 
+    #receiving hosts needs to have open "resource room" for accepting new strains
+    B_PW = b*(
+      h_NU[1] +
+        h_ND[1] + 
+        h_SU[1] +
+        h_SD[1] +
+        h_AU[1] +
+        h_AD[1] 
+    )
+    
+    B_PR = b*(
+      h_NU[2]+
+        h_ND[2]+ 
+        h_SU[2]+
+        h_SD[2]+
+        h_AU[2]+
+        h_AD[2]
+    )
+    
+    
+    
+    # per capita biting rate for hosts who are not treated
+    Bu = (B_PW+B_PR)
+    # per capita biting rate for hosts who are drug treated
+    Bd = B_PR
+    #print(B)
+    
+    # pop size +1 is to prevent N = 0 for undefined scenario
+    fracN = P_acc_N/(NU+ND+1)
+    fracS = P_acc_S/(SU+SD+1)
+    fracA = P_acc_A/(AU+AD+1)
+    
+    # probability that the infected strain has not been seen by the adaptive immune system
+    infN = exp(-fracN/nstrains)
+    infS = exp(-fracS/nstrains)
+    infA = exp(-fracA/nstrains)
+    
+    # Parasite pop ODE------
+    ## sensitive-------
+    dPW_NU.dt = B_PW*NU*infN*(1-rNU/capacity) + #new infections
+      (PW_ND/Tau)- # drug loose effectiveness
+      (PW_NU*Bu*(1-rNU/capacity)*infN*d1) -  #naive class getting infected, and treated
+      (PW_NU*Bu*(1-rNU/capacity)*rho1) - #naive class getting enough # of infections to move to S
+      (PW_NU*Bu*(1-rNU/capacity)*muNU)- #naive class die from infections
+      (PW_NU*att)- #normal death rate
+      (PW_NU*muWU)#clearance of infection
+    
+    
+    #print(dNU.dt)
+    
+    dPW_ND.dt = B_PW*ND*infN*(1-rND/capacity) + #new infections
+      (PW_NU*Bu*(1-rNU/capacity)*infN*d1) - #naive class getting infected, and treated
+      (PW_ND/Tau) - # drug loose effectiveness
+      (PW_ND*Bd*(1-rND/capacity)*rho1) -#naive class getting enough # of infections to move to S
+      (PW_ND*Bd*(1-rND/capacity)*muND) -#naive class die from infections
+      (PW_ND*att)- #normal death rate
+      (PW_ND*muWD)#clearance from drug
+    
+    #print(dND.dt)
+    
+    dPW_SU.dt = B_PW*SU*infS*(1-rSU/capacity) + #new infections
+      (PW_NU*Bu*(1-rNU/capacity)*rho1) +#naive class getting enough # of infections to move to S
+      (PW_SD/Tau)-  # drug loose effectiveness
+      (PW_SU*Bu*(1-rSU/capacity)*infS*d2)- #symptomatic class getting infected, show symptoms, and treated
+      (PW_SU*Bu*(1-rSU/capacity)*rho2) -#symptomatic class getting enough # of infections to move to AS
+      (PW_SU*att)- #normal death rate
+      (PW_SU*muWU)#clearance of infection
+    
+    dPW_SD.dt = B_PW*SD*infS*(1-rSD/capacity) + #new infections
+      (PW_ND*Bd*(1-rND/capacity)*rho1)+ #naive class getting enough # of infections to move to S
+      (PW_SU*Bu*(1-rSU/capacity)*infS*d2)- #symptomatic class getting infected, show symptoms, and treated
+      (PW_SD/Tau) - # drug loose effectiveness
+      (PW_SD*Bd*(1-rSD/capacity)*rho2)- #symptomatic class getting enough # of infections to move to AS
+      (PW_SD*att)- #normal death rate
+      (PW_SD*muWD)#clearance from drug
+    
+    dPW_AU.dt = B_PW*AU*infA*(1-rAU/capacity) + #new infections
+      (PW_SU*Bu*(1-rSU/capacity)*rho2)+ #symptomatic class getting enough # of infections to move to AS
+      (PW_AD/Tau)-  # drug loose effectiveness
+      (PW_AU*Bu*(1-rAU/capacity)*infA*omega*d3)- #asymptomatic class getting infected, show symptoms, and treated
+      (PW_AU*att)- #normal death rate
+      (PW_AU*muWU)#clearance from infection
+    
+    dPW_AD.dt = B_PW*AD*infA*(1-rAD/capacity) + #new infections
+      (PW_SD*Bd*(1-rSD/capacity)*rho2) + #symptomatic class getting enough # of infections to move to AS
+      (PW_AU*Bu*(1-rAU/capacity)*infA*omega*d3)- #asymptomatic class getting infected, show symptoms, and treated
+      (PW_AD/Tau)-  # drug loose effectiveness
+      (PW_AD*att)- #normal death rate
+      (PW_AD*muWD)#clearance from drug
+    
+    ## resistant-------
+    dPR_NU.dt = B_PR*NU*infN*(1-rNU/capacity) + #new infections
+      (PR_ND/Tau)- # drug loose effectiveness
+      (PR_NU*Bu*(1-rNU/capacity)*infN*d1) -  #naive class getting infected, and treated
+      (PR_NU*Bu*(1-rNU/capacity)*rho1) - #naive class getting enough # of infections to move to S
+      (PR_NU*Bu*(1-rNU/capacity)*muNU)- #naive class die from infections
+      (PR_NU*att)- #normal death rate
+      (PR_NU*muRU)#clearance of infection
+    
+    #print(dNU.dt)
+    
+    dPR_ND.dt = B_PR*ND*infN*(1-rND/capacity) + #new infections
+      (PR_NU*Bu*(1-rNU/capacity)*infN*d1) - #naive class getting infected, and treated
+      (PR_ND/Tau) - # drug loose effectiveness
+      (PR_ND*Bd*(1-rND/capacity)*rho1) -#naive class getting enough # of infections to move to S
+      (PR_ND*Bd*(1-rND/capacity)*muND) -#naive class die from infections
+      (PR_ND*att)- #normal death rate
+      (PR_ND*muRD)#clearance from drug
+    
+    #print(dND.dt)
+    
+    dPR_SU.dt = B_PR*SU*infS*(1-rSU/capacity) + #new infections
+      (PR_NU*Bu*(1-rNU/capacity)*rho1) +#naive class getting enough # of infections to move to S
+      (PR_SD/Tau)-  # drug loose effectiveness
+      (PR_SU*Bu*(1-rSU/capacity)*infS*d2)- #symptomatic class getting infected, show symptoms, and treated
+      (PR_SU*Bu*(1-rSU/capacity)*rho2) -#symptomatic class getting enough # of infections to move to AS
+      (PR_SU*att)- #normal death rate
+      (PR_SU*muRU)#clearance of infection
+    
+    dPR_SD.dt = B_PR*SD*infS*(1-rSD/capacity) + #new infections
+      (PR_ND*Bd*(1-rND/capacity)*rho1)+ #naive class getting enough # of infections to move to S
+      (PR_SU*Bu*(1-rSU/capacity)*infS*d2)- #symptomatic class getting infected, show symptoms, and treated
+      (PR_SD/Tau) - # drug loose effectiveness
+      (PR_SD*Bd*(1-rSD/capacity)*rho2)- #symptomatic class getting enough # of infections to move to AS
+      (PR_SD*att)- #normal death rate
+      (PR_SD*muRD)#clearance from drug
+    
+    dPR_AU.dt = B_PR*AU*infA*(1-rAU/capacity) + #new infections
+      (PR_SU*Bu*(1-rSU/capacity)*rho2)+ #symptomatic class getting enough # of infections to move to AS
+      (PR_AD/Tau)-  # drug loose effectiveness
+      (PR_AU*Bu*(1-rAU/capacity)*infA*omega*d3)- #asymptomatic class getting infected, show symptoms, and treated
+      (PR_AU*att)- #normal death rate
+      (PR_AU*muRU)#clearance from infection
+    
+    dPR_AD.dt = B_PR*AD*infA*(1-rAD/capacity) + #new infections
+      (PR_SD*Bd*(1-rSD/capacity)*rho2) + #symptomatic class getting enough # of infections to move to AS
+      (PR_AU*Bu*(1-rAU/capacity)*infA*omega*d3)- #asymptomatic class getting infected, show symptoms, and treated
+      (PR_AD/Tau)-  # drug loose effectiveness
+      (PR_AD*att)- #normal death rate
+      (PR_AD*muRD)#clearance from drug
+    
+    # Generalized immunity ODE------
+    dNU.dt = birth + 
+      (ND/Tau)- # drug loose effectiveness
+      (NU*Bu*(1-rNU/capacity)*infN*d1) -  #naive class getting infected, and treated
+      (NU*Bu*(1-rNU/capacity)*rho1) - #naive class getting enough # of infections to move to S
+      (NU*Bu*(1-rNU/capacity)*muNU)- #naive class die from infections
+      (NU*att)#normal death rate
+    
+    
+    #print(dNU.dt)
+    
+    dND.dt = (NU*Bu*(1-rNU/capacity)*infN*d1) - #naive class getting infected, and treated
+      (ND/Tau) - # drug loose effectiveness
+      (ND*Bd*(1-rND/capacity)*rho1) -#naive class getting enough # of infections to move to S
+      (ND*Bd*(1-rND/capacity)*muND) -#naive class die from infections
+      (ND*att) #normal death rate
+    #print(dND.dt)
+    
+    dSU.dt =  (NU*Bu*(1-rNU/capacity)*rho1) +#naive class getting enough # of infections to move to S
+      (SD/Tau)-  # drug loose effectiveness
+      (SU*Bu*(1-rSU/capacity)*infS*d2)- #symptomatic class getting infected, show symptoms, and treated
+      (SU*Bu*(1-rSU/capacity)*rho2) -#symptomatic class getting enough # of infections to move to AS
+      (SU*att) #normal death rate
+    
+    dSD.dt = (ND*Bd*(1-rND/capacity)*rho1)+ #naive class getting enough # of infections to move to S
+      (SU*Bu*(1-rSU/capacity)*infS*d2)- #symptomatic class getting infected, show symptoms, and treated
+      (SD/Tau) - # drug loose effectiveness
+      (SD*Bd*(1-rSD/capacity)*rho2)- #symptomatic class getting enough # of infections to move to AS
+      (SD*att) #normal death rate
+    
+    dAU.dt = (SU*Bu*(1-rSU/capacity)*rho2)+ #symptomatic class getting enough # of infections to move to AS
+      (AD/Tau)-  # drug loose effectiveness
+      (AU*Bu*(1-rAU/capacity)*infA*omega*d3)- #asymptomatic class getting infected, show symptoms, and treated
+      (AU*att) #normal death rate
+    
+    dAD.dt = (SD*Bd*(1-rSD/capacity)*rho2) + #symptomatic class getting enough # of infections to move to AS
+      (AU*Bu*(1-rAU/capacity)*infA*omega*d3)- #asymptomatic class getting infected, show symptoms, and treated
+      (AD/Tau)-  # drug loose effectiveness
+      (AD*att) #normal death rate
+    
+    #sprintf("%.0f, %.0f,%.0f,%.0f,%.0f,%.0f",dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt)
+    
+    # track the mean number of times infected per immunity class------
+    dP_acc_N.dt = NU*Bu*(1-rNU/capacity) + ND*Bd*(1-rND/capacity) - # total number of infection received
+      (NU*(Bu*(1-rNU/capacity)*muNU+att))*P_acc_N/(NU+ND+1)- #when naive U class die how many infections deducted
+      (ND*(Bd*(1-rND/capacity)*muND+att))*P_acc_N/(NU+ND+1)- #when naive D class die many infections deducted
+      ((NU*Bu*(1-rNU/capacity)*rho1)+(ND*Bd*(1-rND/capacity)*rho1))*P_acc_N/(NU+ND+1) #when naive class leave to S
+    
+    dP_acc_S.dt = ((NU*Bu*(1-rNU/capacity)*rho1)+(ND*Bd*(1-rND/capacity)*rho1))*P_acc_N/(NU+ND+1)+ #when naive class leave to S
+      SU*Bu*(1-rSU/capacity) + SD*Bd*(1-rSD/capacity)- #total number of new infections
+      att*P_acc_S- # normal death
+      ((SD*Bd*(1-rSD/capacity)*rho2)+ (SU*Bu*(1-rSU/capacity)*rho2))*P_acc_S/(SU+SD+1) #Symptomatic move to Asymptomatic
+    
+    dP_acc_A.dt = ((SD*Bd*(1-rSD/capacity)*rho2)+ (SU*Bu*(1-rSU/capacity)*rho2))*P_acc_S/(SU+SD+1)+ #Symptomatic move to Asymptomatic
+      AU*Bu*(1-rAU/capacity) + AD*Bd*(1-rAD/capacity)- #total number of new infections
+      (ImmLost+att)*P_acc_A # normal death
+    
+    a <- 365/(2*pi)
+    db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a))
+    #db.dt =0
+    #sprintf("%.0f, %.0f,%.0f,%.0f",dP_acc_N.dt,dP_acc_S.dt,dP_acc_A.dt,dbeta.dt)
+    
+    return(list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt,
+                  dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt,
+                  dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,
+                  dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt)))
+  })
+}
+
+allPara<-read.csv(paramFile,sep=",",header=T)
+
+currentPara<-allPara[allPara$No==No,]
+
+t2 = seq(from=0,to=100*365,by=1) 
+N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
+       PR_NU=2*(1-WildStrainOnly),PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
+       NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
+       P_acc_N=0,P_acc_S=0,P_acc_A=0,
+       b=currentPara$bavgnull*(1-currentPara$seasonality))
+out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)
+
+out2<-as_tibble(out2)
+sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
+                       PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
+                group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
+
+while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/sumOut2[30,"meanP"]>0.001 )|
+       (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/sumOut2[30,"meanN"]>0.001)){
+  print(out2[100*365,2:23])
+  t2<-t2+30*365
+  N0 = as.numeric(out2[100*365,2:23])
+  names(N0)<-colnames(out2[,2:23])
+  out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)
+  out2<-as_tibble(out2)
+  sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
+                         PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
+    group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
+}
+out2<-as.data.frame(out2)
+out2<-out2%>%mutate(m_N = P_acc_N/(NU+ND),m_S =P_acc_S/(SU+SD), m_A = P_acc_A/(AU+AD),
+                    PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
+                    PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD, 
+                    P=PW+PR,N = NU+ND+SU+SD+AU+AD, prev = 1-exp(-P/N))
+
+##you can check the output through this graph
+#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10000), id.vars=1)
+#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+
+#  facet_wrap(.~variable,scales="free_y")
+
+outFile<-paste(resultsFolder,"/",strsplit(paramFile,split="[.]")[[1]][1],"_results.csv",sep="")
+if(file.exists(outFile)){
+  write.table(as.data.frame(cbind(currentPara,tail(out2,1))),
+            file = outFile, row.names=F, col.names=F,sep=",",
+            quote=F,append=T)
+}else{
+  write.csv(as.data.frame(cbind(currentPara,tail(out2,1))),
+            file = outFile, row.names=F, 
+            quote=F)
+  
+}
+