I am fairly new to python and have a silly question. I looked up for answers, couldnt find it, so here I am (First question on StackOverflow, so forgive if any mistakes).

I have a function (functions.Reservoir) that works perfectly. It takes a complex numpy 2D array and one float value as argument. It returns back another complex numpy 2D array. I try to save all these return values in a new numpy array. The problem is that the first print command gives me the right values, which means my function is working properly. But outside the for loop, when I print the l-variable, all elements have the same value, namely the value of that the last element was supposed to have.

Can you tell me what I'm doing wrong?

another thing you could help me with is a way to save these return values as a 3D numpy array. I don't want to use vstack, because sometimes my loop iterates 4000times and vstack would slow it down massively. What I do now is save the numpy arrays as objects and then later create a new array from this iterable. If there is a better way, please let me know. Thanks in advance

    k= np.zeros((400,3), dtype=np.complex128)
    l= np.empty(noPoints1, dtype=object)
    for i in  xrange(noPoints1):
      l[i] = functions.Reservoir(k, data[ OFFSET1+ Discard +i] )
      print l[i]
      k=l[i]
    print l

Update

Here is more on my code:

def determ(k1, k2,  d): 
    dkdt = np.zeros_like(k1)

    EA1 = np.absolute(k1[0])**2
    EA2 = np.absolute(k1[1])**2

    gain1= g1*(k1[2]-N0)/(1+s*EA1)
    gain2= g2*(k1[2]-N0)/(1+s*EA2)

    dkdt[0]= 0.5 * (1+ A*1j) * (gain1 - pd) * k1[0] 

    dkdt[1]= 0.5 * (1+ A*1j) * (gain2 - pd) * k1[1] 

    dkdt[2]= Stroom/e - k1[2]*ed - EA1*gain1 - EA2*gain2   
    return dkdt

def Reservoir(y,d): 
    q= y.shape[0]               
    mdata= 0.5*np.pi*d*maskvector    #Maskvector 1D with q elements

    for n in xrange(q):  
      yn = y[n,:]               # Predictor

      if n+1<q:
        ytd= y[n+1,:]           # Predictor Delayed Field strength
      else:
        ytd= y[n+1-q,:]
      #See y as the state of a circular reservoir. If you reach the end at q, then it directly connects back at the start.

      fn = determ(k1=yn, k2=ytd, d=mdata[n]) 

      ybar = yn + fn*step           # Corrector
      if n+2<q:
        ytdbar= y[n+2,:]            # Corrector Delayed Field strength  
      else:
        ytdbar= y[n+2-q,:]

      fnbar = determ(k1= ybar, k2=ytdbar, d=mdata[n])

      if n+1<q:
        y[n+1] = yn + 0.5*(fn + fnbar)*step 
      else:
        y[n+1-q] = yn + 0.5*(fn + fnbar)*step 
    return (y)

PS. I tried the answer of HYRY (Collecting results from a loop that returns NumPy Arrays), But still I get the same values for all elements.