Hello,

I have a set of ODEs to model a bacterial infection:

```
function kinInfect!(du, u, p, t)
μ , φ, η, β, ω, χ, ε = p
H = h = 0
#=
du[1] = susceptible
du[2] = infected
du[3] = phages
μ = growth rate (replication coef), a
φ = adsorption rate (transmission coef), b
η = lysis rate, k
β = burst size, L
ω = outflow (decay rate phage), m
H = host response against bacterium
h = host response against phage
=#
x = χ(t, ε)
du[1] = (μ * u[1]) - (φ * u[1] * u[3]) - (H * u[1])
du[2] = (μ * u[2]) + (φ * u[1] * u[3]) - (η * u[2]) - (H * u[2])
du[3] = (x * η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3])
end
Chi(x, y) = (x < y) ? 0.0 : 1.0
# parms
tmax = 20.0 # time span 0-tmax
mu = 0.5 # growth rate (replication coef), a
phi = 1e-7 # adsorption rate (transmission coef), b
eta = 5 # lysis rate, k
beta = 100 # burst size, L
omega = 5 # outflow (decay rate phage), m
tmax = 20.0 # time span 0-tmax
t_in = 2.5 # time of inoculum
v_in = 1e8 # amount of inoculum
s0 = 1000.0 # initial susceptible population
i0 = 0.0 # initial infected population
v0 = 0.0 # initial phage population
tspan = (0.0, tmax)
parms = [mu, phi, eta, beta, omega, Chi, t_in]
u0 = [s0, i0, v0]
# delay infection
condition(u, t, integrator) = t==t_in # time of inoculum
affect!(integrator) = integrator.u[3] += 1e8 # amount of inoculum
cb = DiscreteCallback(condition,affect!)
# implement model
prob = ODEProblem(kinInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])
```

the function Chi is there to neutralize the first term of du[3] until the inoculum time is reached. The first time I ran it I got:

But the second time I got:

I don’t remember having changed anything between the first and second time (apart from removing the vertical line).

How can I debug an ODE function? Are there some techniques?

Thank you