Hello,

I am using DifferentialEquations.jl to solve these ordinary differential equations:

where dN/dt is bacteria susceptible, dI/dt is bacteria infected and dV/dt is phages. I wrote the following code:

```
function dynamo!(du, u, p, t)
μ, κ, φ, ω, η, β = 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
# set parameters
mu = 0.16 # maximum growth rate
kappa = 2.2*10^7 # maximum population density
phi = 10.0^-9 # adsorption rate
beta = 50.0 # burst size
omega = 0.05 # outflow
eta = 1.0 # lyse rate
tmax = 4000.0 # time span 0-tmax
s0 = 50000.0 # initial susceptible population
i0 = 0.0 # initial infected population
v0 = 80.0 # initial phage population
u0 = [s0, i0, v0]
p = [mu, kappa, phi, omega, eta, beta]
tspan = (0.0, tmax)
# set solver
prob = ODEProblem(dynamo!, u0, tspan, p)
soln = solve(prob)
```

and it works, but the plot I get is this:

instead of this:

I guess the difference is in the algorithm chosen by the solver (unless there are flaws in the conversion of the equations).

Can I tune the solver in order to get the expected results? and if I change the algorithm, won’t I introduce a bias in the results, for they will be what they should be but what they are supposed to be? In this case, I have a certain solution available, but what for real life data?

Thank you