I have a population model described by this system of differential equations:

```
dN/dt = μN(1-(N+I)/κ) -φNP -ωN
dI/dt = φNP -ηI -ωI
dP/dt = βηI -φNP -ωP
```

with

```
κ = 2.2 ⋅ 10^7
μ = 0.47
φ = 10^-9
β = 50
η = 1.0
ω = 0.05
```

This system reaches equilibrium after a phase of oscillation. I would like to find the initial value of P that brings N to zero.

To find it, I have equated dN/dt to zero and solved for P, obtaining:

`(μ/φ) * (1-(N+I/κ)) - ω/φ`

N+I is the initial number of susceptible individuals because the infection simply shifts individuals from one compartment to another. I tried it out with this system where I have introduced an extinction clause (if N goes below one, I reset everything):

```
function 𝔅!(du, u, p, t) # SIRl with logistic term
μ, κ, φ, ω, η, β = p
#=
du[1] = susceptible
du[2] = infected
du[3] = phages
=#
du[1] = ((μ * u[1]) * (1 - ((u[1]+u[2])/κ))) - (φ * u[1] * u[3]) - (ω * u[1])
du[2] = (φ * u[1] * u[3]) - (η * u[2]) - (ω * u[2])
du[3] = (β * η * u[2]) - (φ * u[1] * u[3]) - (ω * u[3])
end
# extra parms
tmax = 4000.0 # time span 0-tmax
r_s0 = 50000.0 # initial susceptible population u[1]
r_i0 = 0.0 # initial infected population u[2]
r_v0 = 1000.0 # initial infectious agent population u[3]
Tp = 1000 # time of infection
# execute
Vp = (mu/phi) * (1-(r_s0/kappa)) - omega/phi
tspan = (0.0, tmax)
u0 = [r_s0, r_i0, 0]
parms = [mu, kappa, phi, omega, eta, beta]
# inoculum
condition1(u, t, integrator) = t==Tp
affect1!(integrator) = integrator.u[3] += Vp
cb1 = DiscreteCallback(condition1, affect1!)
# extintion
condition2(u, t, integrator) = u[1]-1
function affect2!(integrator)
integrator.u[1] = 0
integrator.u[2] = 0
end
cb2 = ContinuousCallback(condition2, affect2!)
condition3(u, t, integrator) = u[3]-1
function affect3!(integrator)
integrator.u[3] = 0
end
cb3 = ContinuousCallback(condition3, affect3!)
# run
modification = CallbackSet(cb1, cb2, cb3)
prob = ODEProblem(𝔅!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[Tp])
```

The susceptible N does indeed go down after the infection at t=1000, but does not go extinct:

How can I find the value of P (in terms of concentration and time) that extinguish N? Is there some procedure in DifferentialEquations?

Thank you