timestep variable in ODEINT

Question

First let me apologize if this is a simple question, or it has been answered before. The truth is that I don't even know how to run a search about my question.

Let's say I have the following set of coupled ODE:

dN/dt = (hN(t) + iQ(t))*G(t)

dQ/dt = (jN(t) + lQ(t))*G(t)

This is the code i used to plot the ODEs:

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)

The above does not work, since i assumed i could calculate G(t) for each timestep t, but i noticed len(L) is not the same as Q or N. Basically my question is, how to i calculate a parameter on the ODE that varies with the timestep t?

Thanks a lot