Unexpected solutions solving an ODE using odeint

Question

I am trying to solve a system of 8 coupled differential equations using scipy's odeint. I have already written my code and it runs fine, but the solutions I get are completely different from what I expected. Originally, I would write all these equations within a loop, but since I was having problems with them I wrote them all separately.

The Zd argument that I need to pass is an exponentially decreasing variable the same size as my time vector. However, since odeint only calculates one sample at a time, I need to give it the right Zd sample for each time, so I use a copy of t to get it inside the function. The differential equations have four terms each and they need to be turned on/off at different times, so I use the step function there.

This system of equations is supposed to represent a sort of energy going through adjoining systems, so the eight Z functions are the evolution of energy within each one of these systems over time. At t=0, the energy enters the first of these systems and the others have no energy whatsoever. Zd represents the input energy into the system, which changes over time (it decays exponentially as a result of scattering) and needs to be turned on/off as energy enters/leaves each layer. By the end of my time window, I expected to have little or no energy at all. However, my resulting Zs grow almost exponentially and this doesn't make sense.

I have checked it over and over again but I can't see what is wrong in it. The signs and times in the equations are correct. Any idea on what is wrong here?

This is my code and a plot of the results:

def ODE_solver(Z,t,Q,Zd,Z0,dts,ts,cf,t_copy):

    import numpy as np

    # Define functions to solve for:
    [Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8] = Z
    # Define time intervals:
    [dt1, dt2, dt3, dt4, dt5, dt6, dt7, dt8] = dts    
    # Define cumulative times:
    [t1, t2, t3, t4, t5, t6, t7, t8] = ts    
    # Q factors:
    [Q1, Q2, Q3, Q4, Q5, Q6, Q7, Q8] = Q

    # Set some boundary conditions:
    Z9=0; dt9=1/4 # (It doesn't matter which value I have here, since Z9 = 0)
    t0=0; dt0=1/4; 

    # Find index of the time sample that is being used at the moment:
    t_ind=(np.abs(t-t_copy)).argmin()
    Zd=Zd[t_ind]

    # j=1
    dZ1dt= ((1/(4*dt0)) * Z0 * step_fun (t,t0)) + ((1/(4*dt2)) * Z2 * step_fun (t,t0)) 
    - ((1/(4*dt1)) * Z1 * step_fun (t,t1)) - ((1/(4*dt1)) * Z1 * step_fun(t,t0)) 
    + (((2*np.pi*cf)*Q1) * Zd * step_fun(t,t0) * inv_step_fun(t,t1)) 

    # j=2
    dZ2dt= ((1/(4*dt1)) * Z1 * step_fun (t,t1)) + ((1/(4*dt3)) * Z3 * step_fun (t,t1)) 
    - ((1/(4*dt2)) * Z2 * step_fun (t,t2)) - ((1/(4*dt2)) * Z2 * step_fun(t,t1)) 
    + (((2*np.pi*cf)*Q2) * Zd * step_fun(t,t1) * inv_step_fun(t,t2)) 

    # j=3
    dZ3dt= ((1/(4*dt2)) * Z2 * step_fun (t,t2)) + ((1/(4*dt4)) * Z4 * step_fun (t,t2)) 
    - ((1/(4*dt3)) * Z3 * step_fun (t,t3)) - ((1/(4*dt3)) * Z3 * step_fun(t,t2)) 
    + (((2*np.pi*cf)*Q3) * Zd * step_fun(t,t2) * inv_step_fun(t,t3))     

    # j=4
    dZ4dt= ((1/(4*dt3)) * Z3 * step_fun (t,t3)) + ((1/(4*dt5)) * Z5 * step_fun (t,t3)) 
    - ((1/(4*dt4)) * Z4 * step_fun (t,t4)) - ((1/(4*dt4)) * Z4 * step_fun(t,t3)) 
    + (((2*np.pi*cf)*Q4) * Zd * step_fun(t,t3) * inv_step_fun(t,t4))   

    # j=5
    dZ5dt= ((1/(4*dt4)) * Z4 * step_fun (t,t4)) + ((1/(4*dt6)) * Z6 * step_fun (t,t4)) 
    - ((1/(4*dt5)) * Z5 * step_fun (t,t5)) - ((1/(4*dt5)) * Z5 * step_fun(t,t4)) 
    + (((2*np.pi*cf)*Q5) * Zd * step_fun(t,t4) * inv_step_fun(t,t5))  

    # j=6
    dZ6dt= ((1/(4*dt5)) * Z5 * step_fun (t,t5)) + ((1/(4*dt7)) * Z7 * step_fun (t,t5)) 
    - ((1/(4*dt6)) * Z6 * step_fun (t,t6)) - ((1/(4*dt6)) * Z6 * step_fun(t,t5)) 
    + (((2*np.pi*cf)*Q6) * Zd * step_fun(t,t5) * inv_step_fun(t,t6))  

    # j=7
    dZ7dt= ((1/(4*dt6)) * Z6 * step_fun (t,t6)) + ((1/(4*dt8)) * Z8 * step_fun (t,t6)) 
    - ((1/(4*dt7)) * Z7 * step_fun (t,t7)) - ((1/(4*dt7)) * Z7 * step_fun(t,t6)) 
    + (((2*np.pi*cf)*Q7) * Zd * step_fun(t,t6) * inv_step_fun(t,t7))  

    # j=8
    dZ8dt= ((1/(4*dt7)) * Z7 * step_fun (t,t7)) + ((1/(4*dt9)) * Z9 * step_fun (t,t7)) 
    - ((1/(4*dt8)) * Z8 * step_fun (t,t8)) - ((1/(4*dt8)) * Z8 * step_fun(t,t7)) 
    + (((2*np.pi*cf)*Q8) * Zd * step_fun(t,t7) * inv_step_fun(t,t8))  


    return [dZ1dt,dZ2dt,dZ3dt,dZ4dt,dZ5dt,dZ6dt,dZ7dt,dZ8dt]

def step_fun(x,a):

    if x==a:  return 0.5
    elif x<a: return 0
    elif x>a: return 1 


def inv_step_fun(x,a):

    if x==a:  return 0.5
    elif x<a: return 1
    elif x>a: return 0 

####################################################################################################

import copy
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Number of equations:
N=8

# Define t:
t=np.arange(0,50,0.025)
t_copy=copy.deepcopy(t)
Nt=len(t)

# Define other parameters/arguments:
Q=(0.008017307188069223, 0.008117153569161557, 0.008081825341069762, 0.00803759248462624, 0.00803759248462624, 0.008081825341069762, 0.008117153569161557, 0.008017307188069223)
Z0=np.array([22050, 0., 0., 0., 0., 0., 0., 0.])
Zd0=Z0[0]
dts=np.array([20.67120623,  1.43225437,  1.55738981,  1.63132137,  1.63132137, 1.55738981,  1.43225437, 20.67120623])
ts=np.array([20.67120623, 22.10346059, 23.66085041, 25.29217178, 26.92349315, 28.48088296, 29.91313733, 50.58434356])
cf=1.5

# Define Zd: this is a vector of size Nt, but I will need to give the function only the right sample
# at each time t.
Zd=1/(1+np.exp(t))*5e4

# Solve system of differential equations: 
Z=odeint(ODE_solver,Z0,t,args=(Q,Zd,Zd0,dts,ts,cf,t_copy))

# Let's rearrange this results so it is easier to read in the future:
Zs=[]
for n in range(N):
    Zs.append([])
    for vec in Z:
        Zs[n].append(vec[n])
    Zs[n]=np.array(Zs[n])

plt.figure()
for vec in Zs:
    plt.plot(t,vec)

Answer 1

For the interpretation where the left and right "virtual" compartments are infinite sinks at level Z0=Z9=0, and the walls between the compartments j and j+1 are removed at time ts[j] where also the external filling of compartment j starts, I get the shortened implementation

def ODE_func(Z,t,Q,Zd,Z0,dts,ts,cf,t_copy):

    # Declare some boundary conditions:
    Z9=0; dt9=1/4; t9=ts[-1]+dt9 # (It doesn't matter which value I have here, since Z9 = 0)
    Z0=0; dt0=1/4; t0=0;
    # Define functions to solve for:
    Z = np.concatenate([[Z0],Z,[Z9]]);
    # Define time intervals:
    dt = np.concatenate([[dt0],dts,[dt9]])
    # Define cumulative times:
    ts = np.concatenate([[t0],ts,[t9]]);


    # Find index of the time sample that is being used at the moment:
    Zd=np.interp(t,t_copy,Zd);

    dZdt = [ (Z[j-1]/dt[j-1]-Z[j]/dt[j])*step_fun(t,ts[j-1])/4 
             + (Z[j+1]/dt[j+1]-Z[j]/dt[j])*step_fun(t,ts[j])/4
             + 2*np.pi*cf*Q[j-1]*Zd*step_fun(t,ts[j-1])*step_fun(ts[j],t) for j in range(1,9)]

    return dZdt

and the plots of the "fill heights" Z[j]/dt[j], j=1,...,8 over the long term as

where indeed the levels start to fall after some time.