Переменная временного шага в ODEINT

Question

Сначала позвольте мне извиниться, если это простой вопрос, или на него уже был дан ответ. Правда в том, что я даже не знаю, как провести поиск по моему вопросу.

Допустим, у меня есть следующий набор связанных ODE:

dN / dt = (hN (t) + i Q (t)) * G (t)

dQ / dt = (jN (t) + l Q (t)) * G (t)

Это код, который я использовал для построения ODE:

import numpy 
import matplotlib.pyplot as plt
from scipy import constants
from scipy.integrate import odeint
plt.ion()

#-------------------------------------INITIAL CONSTANT PARAMETERS-------------------------------------------

wavelenght      = 1908e-9 #meter
wp0             = 300e-6 #laser beam diameter
conc            = 1.0 #dopant concentration
crosssection    = 6.5e-14 #spectroscopic absorption cross section in m2
lifetime        = 230e-6 #upper lifetime
cavity_l        = 200e-3 #total lenght of cavity in mm
gain_l          = 10e-3 #lenght of gain medium
n_gain          = 1.82 #refraction index gain medium
R_OC            = 0.3 #reflectivity OC
additional_loss = 0.05 #additional losses
loss_time_max   = 10.0e5 #max loss instroduced by q-switch
r               = 10 #times above threshold laser is pumped before q-switch opens
t0              = 0
tf              = 20e-12
tpulse          = 1e-9
ttotal          = 10e-9

#-------------------------------------INITIAL VARIABLE PARAMETERS-------------------------------------------

#number density from concentration percentage
Ntotalyag             = conc*3*4.55/((3*88.9 + 5*27.0 + 12*16.0)*constants.m_p)
#nominator: 3 at of Y * mass density of Y3Al5O12
#denominator: mass of the Y3Al5O12 unit, calculated from the relative atomic weights and the proton mass

beam_area             = numpy.pi*(wp0/2)*(wp0/2)

roundtrip_time        = 2.0*((cavity_l-gain_l)/constants.c)+(n_gain*gain_l/constants.c) #time for light to travel back and forth inside cavity

losscoef              = - numpy.log(R_OC) + additional_loss 

popinversionthreshold = losscoef * beam_area/(2.0 * crosssection)

Wp = 0.0 #pumping rate

#-------------------------------------FUNCTIONS-------------------------------------------

def f(y, t, L):
    q      = y[0]
    DeltaN = y[1]
    if t > tf:
        losscoef_t = 0.0
    else:
        losscoef_t = loss_time_max - ((loss_time_max - losscoef)/(tf))*t
    # gupta Handbook of photonics ch about micro laser
    L= L.append(losscoef_t)
    f0 = q*((DeltaN/popinversionthreshold)-1.0) * losscoef_t/roundtrip_time
    f1 = Wp*(Ntotalyag-DeltaN) - q*((DeltaN*losscoef_t)/(popinversionthreshold*roundtrip_time)) - DeltaN/lifetime
    return [f0, f1]

#-------------------------------------ODE SOLVING 1-------------------------------------------
# initial conditions
q0 = 0.01               # initial photon population
DeltaN0 = r * popinversionthreshold   # initial inverted population
y0 = [q0, DeltaN0]       # initial condition vector
t  = numpy.linspace(0, ttotal, 1000)   # time grid

L = []

# solve the DEs
soln = odeint(f, y0, t, args=(L,))
Q    = soln[:, 0] 
N    = soln[:, 1] 

#-------------------------------------PLOTS---------------------------------------------------

plt.figure()
plt.plot(t, Q, label='Photon')
plt.plot(t, N, label='Population')
plt.legend(loc=0)

Вышеприведенное не работает, так как я предположил, что могу рассчитать G (t) для каждого временного шага t, но заметил, что len(L) - это не то же самое, что Q или N. В основном мой вопрос состоит в том, как рассчитать параметр для ODE, который меняется с шагом t?

большое спасибо