instead of fixed beta and gama, i want to implement the formula likes this $\backslash$beta $=($Number of New Infections per Unit Time$)/($Number of Susceptible Individuals$)\backslash$times $($Number of Infectious Individuals$),$

$\backslash$gamma $=($Number of Recoveries per Unit Time$)/($Number of Infectious Individuals$),$import numpy as np

import matplotlib.pyplot as plt

# Parameters

beta $=0\mathrm{.}222$ # transmission rate

gamma $=0\mathrm{.}122$ # recovery rate

# Initial conditions

initial$\_$infectious $=$ infectious$\_$data.iloc$[0]$

initial$\_$susceptible $=$ susceptible$\_$data.iloc$[0]$

initial$\_$recoveries $=$ recoveries$\_$data.iloc$[0]$

# Time settings

num$\_$days $=$ len$($dates$)$

h $=10\mathrm{.}5$ # time step

# SIR model differential equations

def sir$\_$model$($t$,$ y$,$ beta, gamma$)$:

S$,$ I, R $=$ y

dSdt $=-$beta $*$ S $*$ I $/$ population

dIdt $=$ beta $*$ S $*$ I $/$ population $-$ gamma $*$ I

dRdt $=$ gamma $*$ I

return $[$dSdt$,$ dIdt, dRdt$]$

# RK$-4$ method

def rk$4\_$method$($t$,$ y$,$ h$,$ beta, gamma$)$:

k$1=$ np$.$array$($sir$\_$model$($t$,$ y$,$ beta, gamma$))$

k$2=$ np$.$array$($sir$\_$model$($t $+$ h$/2,$ y $+$ h$/2*$ k$1,$ beta, gamma$))$

k$3=$ np$.$array$($sir$\_$model$($t $+$ h$/2,$ y $+$ h$/2*$ k$2,$ beta, gamma$))$

k$4=$ np$.$array$($sir$\_$model$($t $+$ h$,$ y $+$ h $*$ k$3,$ beta, gamma$))$

return y $+$ h$/6*($k$1+2*$k$2+2*$k$3+$ k$4)$

# Heun's method

def heun$\_$method$($t$,$ y$,$ h$,$ beta, gamma$)$:

y$\_$predictor $=$ y $+$ h $*$ np$.$array$($sir$\_$model$($t$,$ y$,$ beta, gamma$))$

y$\_$corrector $=$ y $+$ h$/2*($np$.$array$($sir$\_$model$($t$,$ y$,$ beta, gamma$))+$ np$.$array$($sir$\_$model$($t $+$ h$,$ y$\_$predictor, beta, gamma$)))$

return y$\_$corrector

# Initial conditions vector

initial$\_$conditions$\_$rk$4=[$initial$\_$susceptible, initial$\_$infectious, initial$\_$recoveries$]$

initial$\_$conditions$\_$heun $=[$initial$\_$susceptible, initial$\_$infectious, initial$\_$recoveries $+10]$ # Adding $10$ to make initial conditions different

# Initialize arrays to store results

result$\_$rk$4=$ np$.$zeros$(($num$\_$days, $3))$

result$\_$heun $=$ np$.$zeros$(($num$\_$days, $3))$

# Set initial conditions

result$\_$rk$4[0,$ :$]=$ initial$\_$conditions$\_$rk$4$

result$\_$heun$[0,$ :$]=$ initial$\_$conditions$\_$heun

# Time integration loop

for i in range$(1,$ num$\_$days$)$:

t $=$ i $*$ h

result$\_$rk$4[$i$,$ :$]=$ rk$4\_$method$($t$,$ result$\_$rk$4[$i$-1,$ :$],$ h$,$ beta, gamma$)$

result$\_$heun$[$i$,$ :$]=$ heun$\_$method$($t$,$ result$\_$heun$[$i$-1,$ :$],$ h$,$ beta, gamma$)$

print$($f$"$Day $\{$i$\}$: RK$4=\{$result$\_$rk$4[$i$,$ :$]\},$ Heun $=\{$result$\_$heun$[$i$,$ :$]\}")$

# Plot results for RK$-4$

dates $=$ np$.$arange$(1,$ len$($result$\_$rk$4)+1)$ # result$\_$rk$4$ contains data for each day

plt$.$figure$($figsize$=(10,6))$

plt$.$plot$($dates$,$ result$\_$rk$4[$:$,0],$ label$=$'Susceptible $($RK$-4)\text{'})$

plt$.$plot$($dates$,$ result$\_$rk$4[$:$,1],$ label$=$'Infectious $($RK$-4)\text{'})$

plt$.$plot$($dates$,$ result$\_$rk$4[$:$,2],$ label$=$'Recovered $($RK$-4)\text{'})$

# Plot results for Heun

dates $=$ np$.$arange$(1,$ len$($result$\_$heun$)+1)$ #result$\_$heun contains data for each day

plt$.$plot$($dates$,$ result$\_$heun$[$:$,0],$ label$=$'Susceptible $($Heun$)\text{'})$

plt$.$plot$($dates$,$ result$\_$heun$[$:$,1],$ label$=$'Infectious $($Heun$)\text{'})$

plt$.$plot$($dates$,$ result$\_$heun$[$:$,2],$ label$=$'Recovered $($Heun$)\text{'})$

# Customize plot

plt$.$xlabel$(\text{'}$Days$\text{'})$

plt$.$ylabel$(\text{'}$Number of Cases'$)$

plt$.$title$(\text{'}$SIR Model with RK$-4$ and Heun Methods'$)$

plt$.$legend$($loc$=$'upper right'$)$

plt$.$show$()$