Hi all

I have a simple, epidemiological ODE model with three states: S(usceptibles), I(nfected) and R(ecovered). Transmission is modelled with a periodical, seasonal function. Depending on the initial conditions and parameter values, the I state converges towards stable yearly peaks of equal peak size or alternating peak size (see pics).

I would like the integration to stop when this “steady state” is reached, i.e. when the peaks have a consistent size. It isn’t a steady state in the mathematical sense, since the equations depend on time (through the seasonal forcing). I thought about doing with this a continuous callback function to capture the time point when this happens for the first time. The states at this moment would then be used as initial values for another model. However, I have some difficulties with the callback function. The main idea is to monitor the second state (I) and to terminate when its value is equal to its value 2 years (730 days) ago. I realize this is not the best condition to choose, but I couldn’t think of anything better for now.

Below is my attempt, which throws an error. I appreciate any help. Bw, Fabienne

```
using DifferentialEquations, StatsBase, StatsPlots, Random, StatsBase, LabelledArrays, DiffEqCallbacks
Random.seed!(42)
# ODE model: simple SIR model with seasonally forced contact rate
function SIR!(du,u,p,t)
# states
S, I, R = u
N = S + I + R
# params
β = p.β
η = p.η
φ = p.φ
ω = p.ω
μ = p.μ
σ = p.σ
# FOI
βeff = β * (1.0+η*sin(2.0*π*(t-φ)/365.0))
λ = βeff*I/N
# change in states
du[1] = μ*N - λ*S - μ*S + ω*R
du[2] = λ*S - σ*I - μ*I
du[3] = σ*I - μ*R - ω*R
end
# callback function
function condition_terminate(u,t,integrator)
floor(abs(u[2][t]-u[2][t-730]))
end
affect!(integrator) = terminate!(integrator)
cb_terminate = ContinuousCallback(condition_terminate, affect!)
# Solver settings
tspan = (0.0, 7300.0)
abstol = 1.0e-8
reltol = 1.0e-8
maxiters = 1e7
saveat = 1.0
solver = AutoVern7(Rodas5())
# Initiate ODE problem
p = @LArray [0.28,0.07,20,1.0/365,1.0/(80*365),1.0/5.0] (:β,:η,:φ,:ω,:μ,:σ)
u0 = @LArray [9999.0,1.0,0.0] (:S,:I,:R)
# ODE problem and solution
problem = ODEProblem(SIR!,u0,tspan,p)
sol = solve(problem, solver,
abstol=abstol, reltol=reltol,
maxiters=maxiters,
isoutofdomain=(u,p,t)->any(x->x<0.0,u),
callback=cb_terminate,
saveat=saveat)
```

throws the following error:

```
ERROR: LoadError: MethodError: no method matching getindex(::Float64, ::Float64)
Closest candidates are:
getindex(::Real, ::Num) at C:\Users\micky\.julia\packages\Symbolics\h8kPL\src\register.jl:51
getindex(::Real, ::SymbolicUtils.Symbolic) at C:\Users\micky\.julia\packages\Symbolics\h8kPL\src\register.jl:51
getindex(::Number) at number.jl:75
```