Why is odeint time evolving these two differently in python?

Question

I have two programs, one that can take in N coupled ODEs and one that uses 2 coupled ODEs. In the case that I input the 2 same ODEs into both codes, with the same time span, I get different answers. I know the correct answer, so I can deduce that my N many program is wrong.

Here is the code for the 2 equation dedicated one:

# solve the coupled system dy/dt = f(y, t)
def f(y, t):
    """Returns the collections of first-order 
    coupled differential equations"""
    #v11i = y[0]
    #v22i = y[1]
    #v12i = y[2]
    print y[0]

    # the model equations
    f0 = dHs(tRel,vij)[0].subs(v12,y[2])
    f1 = dHs(tRel,vij)[3].subs(v12,y[2])
    f2 = dHs(tRel,vij)[1].expand().subs([(v11,y[0]),(v22,y[1]),(v12,y[2])])

    return [f0, f1, f2]


# Initial conditions for graphing
v110 = 6               
v220 = 6                 
v120 = 4                  
y0 = [v110, v220, v120]       # initial condition vector

sMesh  = np.linspace(0, 1, 10e3)   # time grid

# Solve the DE's
soln = odeint(f, y0, sMesh) 

and here is the N equation dedicated one:

def f(y, t):
    """Returns the derivative of H_s with initial
    values plugged in"""
    # the model equations 
    print y[0]

    for i in range (0,len(dh)):
        for j in range (0,len(y)):
            dh[i] = dh[i].subs(v[j],y[j])

    dhArray = []  
    for i in range(0,len(dh)):
          dhArray.append(dh[i])

    return dhArray

sMesh  = np.linspace(0, 1, 10e3)   # time grid
dh = dHsFunction(t, V_s).expand()
soln = odeint(f, v0, sMesh)

where dHs(tRel,vij) = dHsFunction(t,V_s) i.e. the exact same ODEs. Similarly y0 and v0 are the exact same. But when I print y[0] in the N many case, I get an output of:

6.0
5.99999765602
5.99999531204
5.97655553477
5.95311575749
5.92967598021
5.69527820744
5.46088043467
5.2264826619
2.88250493418
0.53852720647
-1.80545052124
-25.2452277984
-48.6850050755
-72.1247823527
-306.522555124

as opposed to the 2 dedicated case of:

6.0
5.99999765602
5.99999765602
5.99999531205
5.99999531205
5.98848712729
5.98848712125
5.97702879748
5.97702878476
5.96562028875
5.96562027486
5.91961750442
5.91961733611
5.93039037809
5.93039029335
5.89564277275
5.89564273736
5.86137647436
5.86137638807
5.82758984835

etc.

where the second result is the correct one and graphs the proper graphs.

Please let me know if more code is needed or anything else. Thanks.

Answer 1

Your second version for f modifies the value of the global variable dh. On the first call, you substitute in values in it, and these same values are then used in all subsequent calls.

Avoid that by using e.g. dh_tmp = list(dh) inside the function.