在每次迭代期间，您需要将“a”、“B”、“s”和“w”值存储在单独的列表中。

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

# parameters

S = 0.0001
M = 30.03
K = 113.6561
Vr = 58
R = 8.3145
T = 298.15
Q = 0.000133
Vp = 0.000022
Mr = 36
Pvap = 1400
wf = 0.001
tr = 1200
mass = 40000

# define t
time = 14400
t = np.arange(0, time + 1, 1)

# define initial state
Cv0 = (mass / Vp) * wf  # Cv(0)
Cr0 = (mass / Vp) * (1 - wf)
Cair0 = 0  # Cair(0)

# define function and solve ode
def model(x, t):
    C = x[0]  # C is Cair(t)
    c = x[1]  # c is Cv(t)
    a = Q + (K * S / Vr)
    b = (K * S * M) / (Vr * R * T)
    s = (K * S * M) / (Vp * R * T)
    w = (1 - wf) * 1000
    Peq = (c * Pvap) / (c + w * c * M / Mr)
    Pair = (C * R * T) / M
    dcdt = -s * (Peq - Pair)
    if t <= tr:
        dCdt = -a * C + b * Peq
    else:
        dCdt = -a * C
    return [dCdt, dcdt]

# Initialize lists for a, b, s, and w values
a_values = []
b_values = []
s_values = []
w_values = []

# Solve the ODE for each time step and store a, b, s, and w values
x = np.zeros((len(t), 2))
x[0] = [Cair0, Cv0]

for i in range(1, len(t)):
    x[i] = odeint(model, x[i - 1], [t[i - 1], t[i]])[-1]

    a = Q + (K * S / Vr)
    b = (K * S * M) / (Vr * R * T)
    s = (K * S * M) / (Vp * R * T)
    w = (1 - wf) * 1000

    a_values.append(a)
    b_values.append(b)
    s_values.append(s)
    w_values.append(w)

# Extract the results
C = x[:, 0]
c = x[:, 1]