title: How do I pass function as an array element
body: <p>Here is the updated post: I am providing the entire code. In this code, I is constant. I am trying to make I as a Gaussian function that can be generated using this code:</p>
  <p>Code for Gaussian:</p>
  <pre><code>mean = 0
  
  standard_deviation = 1
  x_values = np.arange(-5, 5, 0.1)
  y_values = scipy.stats.norm(mean, standard_deviation)
  plt.plot(x_values, y_values.pdf(x_values))
  </code></pre>
  <p>Main Code where I want to change constant I to the curve above i.e. I is not a constant anymore but a Gaussian curve in the form I(t):</p>
  <pre><code>from scipy import *
  
  from scipy.integrate import ode
  
  from pylab import *
  
  import matplotlib.pyplot as plt
  
   
  
  ### Simualtion Outputs ###
  
  N = []      # y[0] Carrier concentration
  
  S = []      # y[1] Photon concentration
  
   
  
   
  
   
  
  ### Simualtion input  parameters ###
  
  I = 2.0e-2                                      # p[0] Pumping current (A)
  
  q = 1.6e-19                                     # p[1] Electron charge (C)
  
  V = 9e-11                                       # p[2] Device volume (cm^3)
  
  tn = 1.0e-9                                     # p[3] Carrier relaxation time in seconds (s)
  
  g0 = 1.5e-5                                     # p[4] Gain slope constant (cm^3/s)
  
  Nth = 1e18                                      # p[5] Threshold carrier density (cm^-3)
  
  EPS = 1.5e-17                                   # p[6] Gain compression factor (cm^3)
  
  Gamma = 0.3                                     # p[7] Confinement factor
  
  tp = 1.0e-12                                    # p[8] Photon lifetime in cavity (s)
  
  Beta = 1.0e-4                                   # p[9] Spontaneous Emission Factor
  
   
  
  ### Define equations to be solved
  
  def Laser_rates(t, y, p):
  
     
  
      I = p[0]; q = p[1]; V = p[2]; tn = p[3]; g0 = p[4]; Nth = p[5]; EPS = p[6];  Gamma = p[7]; tp = p[8]; Beta = p[9]
  
      dy = zeros([2])
  
      dy[0] = (I/(q* V)) - (y[0]/tn) -  g0*(y[0] - Nth)*(y[1]/(1 + EPS* y[1]))
  
      dy[1] = Gamma* g0* (y[0] - Nth)*(y[1]/(1 + EPS* y[1])) - y[1]/tp + (Gamma* Beta* y[0]) / tn
  
     
  
      return dy
  
     
  
   
  
  ### Time and initial conditions ### 
  
  t0 = 0; tEnd = 1e-9; dt = 1e-12             #Time constraints
  
  y0 = [4.78e18,1e12]                                 #Initial conditions
  
  Y=[]; T=[]                                  #Create empty lists
  
   
  
  p = [41.8742e-2, 1.6e-19,9e-11, 5e-9,  1.5e-5, 5.75e18,   1.5e-17,  0.2,     1.33e-12,   0.0001]             #Parameters for odes
  
  #p = [ I,    q,     V,   tn,    g0,      Nth,    EPS,    Gamma,     tp,    Beta]
  
   
  
  #Setup integrator with desired parameters
  
  r = ode(Laser_rates).set_integrator('vode', method = 'bdf')
  
  r = ode(Laser_rates).set_integrator('dopri5', nsteps = 1e4)
  
  r.set_f_params(p).set_initial_value(y0, t0)
  
   
  
  ### Simualtion check ###
  
  while r.successful() and r.t+dt &lt; tEnd:
  
      r.integrate(r.t + dt)
  
      Y.append(r.y)        #Makes a list of 1d arrays
  
      T.append(r.t)
  
   
  
  ### Format output ###
  
  Y = array(Y)             #Convert from list to 2d array
  
  N = Y[:, 0]
  
  S = Y[:, 1]
  
   
  
   
  
  ### Plotting ###
  
  f, axarr = plt.subplots(2, sharex=True) # Two subplots, the axes array is 1-d
  
  axarr[0].plot(T, N)
  
  axarr[0].set_ylabel(&quot;Carrier Conc (cm^-3)&quot;)
  
  axarr[0].set_title('Laser-Rate Simulation')
  
  axarr[1].plot(T, S)
  
  axarr[1].set_ylabel(&quot;Photon Conc (cm^-3)&quot;)
  
  axarr[1].set_xlabel(&quot;Time (s)&quot;)
  
  show()
  
   
  </code></pre>