I am trying to do the bifurcation diagram for a nonautonomous SIR two-patch model. But It’s not giving the proper result. The model and parameter values are in the code. I solved the model by ODE solution for a set of bifurcation parameter values and then took the maximum and minimum for the plot. How should I draw a phase plane analysis for my model? Any suggestions? 

Bifurcation Analysis #########

library(deSolve)
library(tidyverse)
vals_betaAbar <- seq(1,70,by=1)  #bifurcation parameter

outdf = NULL

for (i in 1:length(vals_betaAbar)) {

betaAbar = vals_betaAbar[i]


sir.two.patch.model.closed <- function (t, x, params){
#create local variable IA, IB, RA, RB, SA, SB
IA <- x[1]
IB <- x[2]
RA <- x[3]
RB <- x[4]
SA <- x[5]
SB <- x[6]
#we can simplify code using “with”
#this argument to “with” lets us use the variable Names
with(
as.list(params),
{ #the system of rate equations
dIA <- betaAbar(1+sigAcos(2pit/Ti))IASA/Na-(muA+gamaA+phiIAB(1+sigAcos(2pit/Ti)))IA+phiIBA(1+sigAcos(2pit/Ti))IB
dIB <- betaBbar(1+sigAcos(2pit/Ti))IBSB/NB-(muB+gamaB+phiIBA(1+sigAcos(2pit/Ti)))IB +phiIAB(1+sigAcos(2pit/Ti))IA
dRA <- gamaAIA-(muA+phiRAB(1+sigAcos(2pit/Ti)))RA+phiRBA(1+sigAcos(2pit/Ti))RB
dRB <- gamaBIB-(muB+phiRBA(1+sigAcos(2pit/Ti)))RB+phiRAB(1+sigAcos(2pit/Ti))RA
dSA <- muANa-betaAbar(1+sigAcos(2pit/Ti))IASA/Na-(muA+phiSAB(1+sigAcos(2pit/Ti)))SA+phiSBA(1+sigAcos(2pit/Ti))SB
dSB <- muBNB-betaBbar(1+sigAcos(2pit/Ti))IBSB/NB-(muB+phiSBA(1+sigAcos(2pit/Ti)))SB+phiSAB(1+sigAcos(2pit/Ti))SA
#combine results into a single vector dx
dx <- c(dIA,dIB,dRA,dRB,dSA,dSB)
#return result as a list
list(dx)
}
)
}
#function seq returns a sequence
times <- seq(0,100,length.out=3000)
#function c “c”combines values into a vector
params <- c(Na=100, NB=1000, betaAbar=betaAbar, betaBbar=5,muA=0.06, muB=0.02, gamaA=12, gamaB=12,phiIAB=0.4, phiIBA=0.04, phiSAB=0.4, phiSBA=0.04, phiRAB=0.4, phiRBA=0.04, sigA=0., Ti=1)

#initial conditions
xstart <- c(IA=5,IB=5,RA=0,RB=0,SA=100-5,SB=1000-5)
#result stored in data frame
out <- as.data.frame(ode(xstart,times,sir.two.patch.model.closed,params))

out$betaAbar= vals_betaAbar[i]

out=out[251:300,]

outdf = rbind.data.frame(outdf, out)

}


bifdata=NULL
for (i in 1:length(vals_betaAbar)) {

betaAbar = vals_betaAbar

bifdata_A = NULL
bifdata_A$IA_max=max(outdf[outdf$betaAbar==vals_betaAbar[i],]$IA)
bifdata_A$IA_min=min(outdf[outdf$betaAbar==vals_betaAbar[i],]$IA)
bifdata_A$IB_max=max(outdf[outdf$betaAbar==vals_betaAbar[i],]$IB)
bifdata_A$IB_min=min(outdf[outdf$betaAbar==vals_betaAbar[i],]$IB)

bifdata_A$SA_max=max(outdf[outdf$betaAbar==vals_betaAbar[i],]$SA)
bifdata_A$SA_min=min(outdf[outdf$betaAbar==vals_betaAbar[i],]$SA)
bifdata_A$SB_max=max(outdf[outdf$betaAbar==vals_betaAbar[i],]$SB)
bifdata_A$SB_min=min(outdf[outdf$betaAbar==vals_betaAbar[i],]$SB)

bifdata_A$RA_max=max(outdf[outdf$betaAbar==vals_betaAbar[i],]$RA)
bifdata_A$RA_min=min(outdf[outdf$betaAbar==vals_betaAbar[i],]$RA)
bifdata_A$RB_max=max(outdf[outdf$betaAbar==vals_betaAbar[i],]$RB)
bifdata_A$RB_min=min(outdf[outdf$betaAbar==vals_betaAbar[i],]$RB)

bifdata_A$betaAbar= vals_betaAbar[i]

bifdata = rbind.data.frame(bifdata, bifdata_A)
}

dev.new()
ggplot(data=bifdata)+
geom_line(aes(x=betaAbar,y=IA_max))+
geom_point(aes(x=betaAbar,y=IA_min))

#############################