Epidemiological model

I am trying to do this problem (it is part of a bigger set of problems):
The β parameter describes the probability that a contagious person (in E2, I, Ia) meets and infects a susceptible person. In reality, β depends on numerous factors, including the infectiousness of the disease itself, immunity of the population resulting from vaccination and
previous infections, as well as the general behavior of the population. We will now extend our model to use a piecewise constant β. Epidemiologists often refer to the reproduction number R of an epidemic, which is the average number of new persons that an infected person infects. The critical number is R = 1, since if R < 1 the epidemic will decline, while for R > 1
it will grow exponentially. In the simplest models, the relationship between R
and β is R = βτ , where τ is the mean duration of the infectious period. In our
model, which has multiple infectious categories, we have
R = re2β/λ2 + riaβ/µ + β/µ, since the mean durations of the E2 period is 1/λ2 and the mean duration of both I and Ia is 1/µ. The choice of β = 0.4 gives R ≈ 3.2, which is similar to
the values used by the Institute of Public Health (FHI) to model the early stage
of the outbreak in Norway, from mid February to mid March 2020. As we all
know, severe restrictions on travel and social interactions were introduced in
Norway on March 12, 2020. These restrictions substantially reduced the number
of contacts between infected and susceptible persons, which is represented in the
model as a reduction in the parameter β. Write a function beta(t) which represents the piecewise constant β:
β(t) = 
0.4 for t < 30
0.083 for t > 30
Call the function solve_SEIR from Problem 2 with the new piecewise constant
β.

This is my attempt:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from SEIR import SEIR

def beta(t):
    if t < 30:
        return 0.4
    else:
        return 0.083

def solve_SEIR(T, S_0, E2_0, beta):
    model = SEIR(beta=beta)
    U0 = [S_0, 0, E2_0, 0, 0, 0]
    t_eval = np.arange(0, T + 1)
    sol = solve_ivp(model, [0, T], U0, t_eval=t_eval, method='RK45')
    return sol.t, sol.y

def plot_SEIR(t, u, components=['S', 'I', 'Ia', 'R']):
    labels = ['S', 'E1', 'E2', 'I', 'Ia', 'R']
    indices = [labels.index(comp) for comp in components]
    for i in indices:
        plt.plot(t, u[i], label=labels[i])
    plt.legend()
    plt.title('Disease dynamics predicted by the SEIR model')
    plt.xlabel('Time (days)')
    plt.ylabel('Population')
    plt.show()

T = 150  # Total time duration
S_0 = 5.5e6  # Initial susceptible population (approximate population of Norway)
E2_0 = 100  # Initial exposed population in the second latent compartment

t, u = solve_SEIR(T, S_0, E2_0, beta)

plot_SEIR(t, u)

However my plot just displays two flat lines. What is wrong?

Sorry, but could you re-format your code to use the triple-backquotes to start and end the code section, and also fix the indentation? We don’t know how your code actually works without this, so we can’t really help you with the information you’ve given us.

Thanks for the notice! I think it is fixed now:)