I am working on a mathematical model that describes the growth of 4 different bacterial populations and cells of the immune system under certain conditions. The model is governed by the equations below.

**POPULATION 1:**

$ { underbrace { frac {dN_ {PS}} {dt}} _ { text {Rate of variation of population}} = underbrace {r N_ {PS}} {{exponential growth}} cdot underbrace { left (1- frac {N_ {PS} + N_ {PR}} {K} right)} _ { text {Growth Limitation}} – underbrace { theta_ {PS} N_ {PS}} { text {Natural Death}} – underbrace {A_ {PS} N_ {PS}} {{text} Biofilm Formation}} + underbrace {D_ {BS} N_ {BS}} _ { text {Biofilm Dispersion}} – under {} {PS} N_ {PS}} { text {Mutation Rate}} – underbrace { eta delta_ {PS} A_ {m} N_ { PS}} _ { text {Antibiotic Inhibition}} – left { underbrace { Gamma N_ {PS} I} _ { text {Immune System}} right }} $

**POPULATION 2:**

$ { frac {dN_ {BS}} {dt} = (r-c_ {b}) N {{BS} cdot left (1- frac {N_ {BS} + N_ {BR}} {K} right) – theta_ {BS} N_ {BS} + A_ {PS} N_ {PS} -D_ {BS} N_ {BS} – phi_ {BS} N_ {BS} – anda delta_ {BS} A_ {m} N_ {BS} – left { Gamma N_ {BS} I right }} $

**POPULATION 3:**

$ { frac {dN_ {PR}} {dt} = (r-c_ {r}) N {{PR} cdot left (1- frac {N_ {PS} + N_ {PR}} {K} right) – theta_ {PR} N {{PR} -A_ {PR} N {{PR} + D_ {BR} N {{BR} + phi_ {PS} N {{PS} – eta delta_ {PR} A_ {m} N_ {PR} – left { Gamma N_ {PR} I right }} $

**POPULATION 4:**

$ { frac {dN_ {BR}} {dt} = (r- (c_ {b} + c_ {r})) N_ {BR} cdot left (1- frac {N_ {BS} + N_ { BR}} {K} right) – theta_ {BR} N {{}} + A_ {PR} N {{PR} -D_ {BR} N {{}} + phi_ {BS} N_ {BS} – eta delta_ {BR} A_ {m} N_ {BR} – left { Gamma N_ {BR} I right }} $

**NAIV CELLS (IMMUNE SYSTEM):**

$ { frac {dV} {dt} = frac {- sigma VB} { pi + B}} $

**EFFECTOR CELLS (IMMUNE SYSTEM):**

$ {{ frac {dE} {dt} = (2 sigma V + sigma E) cdot frac {B} { pi + B} -hE left (1- frac {B} { pi + B} right)}} $

**MEMORY CELLS (IMMUNE SYSTEM):**

$ {{ frac {dM} {dt} = fEh left (1- frac {B} { pi + B} right)}} $

**TOTAL BACTERIAL DENSITY:**

$ { frac {dB} {dt} = N_ {PS} + N_ {PR} + N_ {BS} + N_ {BR}} $

**DENSITY OF THE IMMUNE SYSTEM:**

$ { frac {dI} {dt} = V + E + M} $

**ANTIBIOTIC UPTAKE 1:**

$ {{ eta} = begin {cases}

1 & t_ {1} leq t leq t_ {1} + t_ {2} \

0 & t_ {1}> t : or : t> t_ {1} + t_ {2}

end {cases}} $

**ANTIBIOTIC UPTAKE 2:**

$ {{ eta} = begin {cases}

1 & B geq varOmega \

0 & B < varOmega

end {cases}} $

I am interested in how $ N_ {PS}, N_ {PR}, N_ {BS}, N_ {BR}, V, E, M $ change over time for which I implemented an algorithm in Python to solve the equations. Greek letters and other undescribed parameters (eg, r, K, etc.) are mainly constants defined at the beginning of the program execution.

An example of the functions used is shown below. Currently, as you can also see in the code, I use an Euler method to solve the equations. However, I would like to implement at least Heun's method or even a higher order Runge-Kutta method to solve them.

I'm already stuck with Heun's method and I do not know how to implement it. I ask for help on how to modify the following code and replace it for example by Heun if that is possible with these equations.

```
def sensitive_bacteria_PS (previous_density, growth_rate, density_PR, maximum_density_P, death_rate_PS, attachment_rate_PS, mutation_rate_PS, antibiose_pol,
previous_density + ((previous_density * (1) (previous_density + density_PR) / maximum_density_P)
```