SEIR

=SEIR(s_initial, e_initial, i_initial, r_initial, beta, gamma, sigma, t_start, t_end, [timesteps], [solve_ivp_method])

=SEIR(990.0, 0.0, 10.0, 0.0, 0.3, 0.1, 0.2, 0.0, 10.0, 10)

=SEIR(990.0, 0.0, 10.0, 0.0, 0.3, 0.1, 0.2, 0.0, 10.0, 10, "RK23")

=SEIR(990.0, 0.0, 10.0, 0.0, 0.3, 0.1, 0.2, 0.0, 50.0, 10)

=SEIR(950.0, 20.0, 30.0, 0.0, 0.5, 0.2, 0.3, 0.0, 10.0, 10)

from scipy.integrate import solve_ivp
import numpy as np
 
def seir(s_initial, e_initial, i_initial, r_initial, beta, gamma, sigma, t_start, t_end, timesteps=10, solve_ivp_method='RK45'):
    """
    Numerically solves the SEIR system of ordinary differential equations for infectious disease modeling.
 
    Args:
        s_initial: Initial number of susceptible individuals.
        e_initial: Initial number of exposed individuals.
        i_initial: Initial number of infected individuals.
        r_initial: Initial number of recovered individuals.
        beta: Infection rate parameter.
        gamma: Recovery rate parameter.
        sigma: Rate at which exposed individuals become infectious.
        t_start: Start time of integration.
        t_end: End time of integration.
        timesteps: Number of timesteps to solve for. Default is 10.
        solve_ivp_method: Integration method ('RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'). Default is 'RK45'.
 
    Returns:
        2D list with header row: t, S, E, I, R. Each row contains time and compartment values, or an error message (str) if input is invalid.
 
    This example function is provided as-is without any representation of accuracy.
    """
    # Validate inputs
    try:
        s0 = float(s_initial)
        e0 = float(e_initial)
        i0 = float(i_initial)
        r0 = float(r_initial)
        beta_ = float(beta)
        gamma_ = float(gamma)
        sigma_ = float(sigma)
        t0 = float(t_start)
        t1 = float(t_end)
        ntp = int(timesteps)
    except Exception:
        return "Invalid input: all initial values, parameters, and timesteps must be numbers."
    if t1 <= t0:
        return "Invalid input: t_end must be greater than t_start."
    if ntp <= 0:
        return "Invalid input: timesteps must be a positive integer."
    if any(x < 0 for x in [s0, e0, i0, r0, beta_, gamma_, sigma_]):
        return "Invalid input: all initial values and parameters must be non-negative."
    allowed_methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']
    if solve_ivp_method not in allowed_methods:
        return f"Invalid input: solve_ivp_method must be one of {allowed_methods}."
    # Create time array for evaluation
    t_eval = np.linspace(t0, t1, ntp)
    
    N = s0 + e0 + i0 + r0
    def seir_ode(t, y):
        S, E, I, R = y
        dSdt = -beta_ * S * I / N
        dEdt = beta_ * S * I / N - sigma_ * E
        dIdt = sigma_ * E - gamma_ * I
        dRdt = gamma_ * I
        return [dSdt, dEdt, dIdt, dRdt]
    try:
        sol = solve_ivp(
            seir_ode,
            [t0, t1],
            [s0, e0, i0, r0],
            method=solve_ivp_method,
            dense_output=False,
            rtol=1e-6,
            atol=1e-8,
            t_eval=t_eval
        )
    except Exception as e:
        return f"ODE integration error: {e}"
    if not sol.success:
        return f"ODE integration failed: {sol.message}"
    result = [["t", "S", "E", "I", "R"]]
    for i in range(len(sol.t)):
        t_val = float(sol.t[i])
        s_val = float(sol.y[0][i])
        e_val = float(sol.y[1][i])
        i_val = float(sol.y[2][i])
        r_val = float(sol.y[3][i])
        # Disallow nan/inf
        if any([
            not isinstance(t_val, float) or not isinstance(s_val, float) or not isinstance(e_val, float) or not isinstance(i_val, float) or not isinstance(r_val, float),
            t_val != t_val or s_val != s_val or e_val != e_val or i_val != i_val or r_val != r_val,
            abs(t_val) == float('inf') or abs(s_val) == float('inf') or abs(e_val) == float('inf') or abs(i_val) == float('inf') or abs(r_val) == float('inf')
        ]):
            return "Invalid output: nan or inf detected."
        result.append([t_val, s_val, e_val, i_val, r_val])
    return result