I have a set of ODEs to model a bacterial infection:
function kinInfect!(du, u, p, t) μ , φ, η, β, ω, χ, ε = p H = h = 0 #= du = susceptible du = infected du = 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 = (μ * u) - (φ * u * u) - (H * u) du = (μ * u) + (φ * u * u) - (η * u) - (H * u) du = (x * η * β * u) - (φ * u * u) - (ω * u) - (h * u) 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 += 1e8 # amount of inoculum cb = DiscreteCallback(condition,affect!) # implement model prob = ODEProblem(kinInfect!, u0, tspan, parms) soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=)
the function Chi is there to neutralize the first term of du 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?