The resolution was coded in python language using the scipy package. The function finally applied was solve_ivp under LSODA condition. This is a function already implemented in python in order to numerically integrate a system of ordinary differential equations given initial values.
We studied the equilibrium state reach for each type of infection depending on the load capacity, the initial number of susceptible bacteria and of phages injected. The goal is also to see the metabolic impact of the phage infection on the bacteria.

Table 1: Graphs obtained depending on the values of N*, S, φ, and A for chronic infection
Type of infection N* S φ A Graph
chronic 2e6 1e6 1e4 1000
chronic 2e6 1e5 1e5 1000
chronic 2e4 1e4 1e3 1000

These three simulations show the impact of the initial susceptible population and N*. The metabolic cost is not studied as A is fixed at 1000 for each run.
The two parameters determine the speed at which the population is infected. In this chronic model the phage always takes control of the bacteria population that reach the load capacity.
The results are not shown here but the value of A does not invert the population and the speed is almost not impacted by a high value of A.

Table 2: Graphs obtained depending on the values of N*, S, φ, and A for lytic infection
Type of infection N* S φ A Graph
lytic 2e6 1e6 1 1000 and 0 (same graph)
lytic 2e5 1e5 1 1000
lytic 2e5 1e5 1 0
lytic 2e4 1e4 100 1000 and 0

This time we study the lytic infection, there are two types of equilibrium depending of the initial values of N*and S.
First we study oscillations with less and less bacteria and phages, this result is obtained whatever the initial number of phage is for these values of A and S. The second equilibrium is a stabilization of the populations, the susceptible population surpassing the infected one.
Here the metabolic cost has an impact on the bacteria growth but not on the type of equilibrium.