# How to set a systems of ODEs to zero with Boundary Value Problem?

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 = susceptible
=#
du = ((μ * u) * (1 - ((u+u)/κ))) - (φ * u * u) - (ω * u)
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 = susceptible
du = infected
du = phages
=#
du = ((μ * u) * (1 - ((u+u)/κ))) - (φ * u * u) - (ω * u)
du = (φ * u * u) - (η * u) - (ω * u)
du = (β * η * u) - (φ * u * u) - (ω * u)
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

I tried with:

``````julia> function bc!(residual, u, p, t)
residual = (u - N0)/N0 # keep the starting value constant
residual = (u - I0)/I0 # keep the starting value constant
residual = u[end]/N0 # u[end] should be zero at the final timestep
end
bc! (generic function with 1 method)

julia> N0 = 1.0*10^7
1.0e7

julia> I0 = 1E5
100000.0

julia> P0 = 2*10^8
200000000

julia> pGuess = 1.0
1.0

julia> tspan = (0.0,1500.0)
(0.0, 1500.0)

julia> bvp0 = BVProblem(modelFun!, bc!, [N0,I0,P0], tspan)
BVProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 1500.0)
u0: [1.0e7, 100000.0, 2.0e8]

julia> sol0 = solve(bvp0, Shooting(Vern7()), dtmin=0.01, dt=0.5)
ERROR: MethodError: no method matching iterate(::DiffEqBase.NullParameters)
Closest candidates are:
iterate(::Base.EnvDict) at env.jl:119
iterate(::Base.EnvDict, ::Any) at env.jl:119
iterate(::JuliaInterpreter.ExprSplitter) at /home/gigiux/.julia/packages/JuliaInterpreter/LBEYA/src/construct.jl:535
...
``````