I have a system of ordinary differential equations (ODEs) that model the interaction between a bacterium and phage:

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

I am looking for the values of P that bring N to zero in a stable term. As is, the addition of P to the system reduces N, but then N and P go into oscillation and then they reach an equilibrium where N >> 0.

I understand this might be a problem for Boundary Value Problem but how do I set the system?

I understand BVP requires a function to work upon, which probably is the dN/dt I indicated, hence:

```
function fun!(du, u, p, t)
μ, κ, φ, ω = p
#=
du[1] = susceptible
=#
du[1] = ((μ * u[1]) * (1 - ((u[1]+u[2])/κ))) - (φ * u[1] * u[3]) - (ω * u[1])
end
```

(do I need to include dI/dt?) This is derived from the system of ODEs I set for the model:

```
function modelFun!(du, u, p, t) # base infection model 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
```

I then need an implicit function bc, but what would that be?

How is the setting of this kind of problem?

The parameters I used are:

```
κ = 2.2 ⋅ 10^7
μ = 0.47
φ = 10^-9
β = 50
η = 1.0
ω = 0.05
N0 = 1.0*10^7
P0 = 2*10^8
```

Thank you