How to solve two first order ODEs in which one of the variables contains a list of elements?

Question

I need help with solving three first-order ODEs using the scipy.integrate.ode module with the integration method of the Runge Kutta Method. My problem is that I don't know how to solve the ODE with one of the variables as a list of 1000 elements. As shown in the code section below, I want to loop through this list of 1000 values one by one when each time the ODEs are being integrated. The values "iI, e, V, Ts, Vg, Ag, N0, Tp, Gamma, Device_length" are all fixed constants. The value of "Pin_dbm" should change when each time the ODE solver is called. The unknown variables that are needed to be solved are "y[0]", "y[1]", and "dy[2]". How do I do this?

This is what I have tried:

    # Import libraries
    import matplotlib.pyplot as plt
    import numpy as np
    import os
    from scipy.integrate import ode

    # User input
    I = int(input('What is the SOA pumping current [mA]: '))        
    Device_length = float(input('What is the device (SOA) length [m]: '))

    # Fixed value constants
    iI = I*1e-3             
    Ts = 1e-9               
    Tp = 1.5e-12            
    Vg = 9e7                
    e = 1.6e-19             
    V = 50e-4*1.8e-15      
    Gamma = 0.052           
    Ag = 3e-20              
    N0 = 5e18        
    Mode  = 1       
    
    # 1000 elements list
    Pin = np.arange(1,100,0.1, dtype=float) 
    Pin_dbm = [10*np.log(P) for P in Pin]    # the list of 1000 elements
    L = np.arange(0,50,0.5, dtype=float)

    N = []                  
    S = []                  
    G = []                  
    N_end = []               
    S_end = []              
    DL = []

    def solver(Pin_dbm):
        global N
        global S
        global G

    # Function definition
        def RK4_ODE_Solver(t, y, p):

            dy = np.zeros([3])

            # three ODEs
            dy[0] = (iI/e*V) - (y[0]/Ts) - Vg*Ag*(y[0]-N0)*y[1]               
            dy[1] = y[1]/Tp + Vg*Gamma*Ag*(y[0] - N0)*y[1] + Pin_dbm          
            dy[2] = 10*np.log(np.exp(Gamma*Ag*(y[0] - N0)*Device_length))

            return dy
    
            # Initial time conditions
            t0 = 0; t_end = 1e-8; dt = 1e-13               
            y_initial = [1e16, 0]                         
            Y = [];                                
            p = [iI, e, V, Ts, Vg, Ag, N0, Tp, Gamma, [i for i in Pin_dbm]]       

            # Using RUNGE KUTTA 4th order to integrate the ODEs
            r = ode(RK4_ODE_Solver).set_integrator('dopri5', nsteps = 1e-4)
            # r = RK45(SOA_rates, t0, y_initial, t_bound=t_end-t0, first_step=1e-4)
            # r = solve_ivp(SOA_rates, method='RK45', t_span=(t0, t_end), y0=y_initial, dense_output=True) 
            r.set_f_params(p).set_initial_value(y_initial, t0)

            while r.successful() and r.t + dt < t_end:
                r.integrate(r.t + dt)
                Y.append(r.y)

            Y = np.array(Y)
            N = Y[:, 0]
            S = Y[:, 1]
            G = Y[:, 2]

            S_end.append(S[-1:])
            N_end.append(N[-1:])

    # Main function
    if (Mode == 1):
        Solver(Pin_dbm)