Use numpy arrays as arguments in odeint

I am trying to solve a system with differential equations using odeint. I have 4 txt files (that look like the picture below). I read them and I save them in numpy arrays (length:8000) (maby not with the most effective way, but anyway...). I want to pass these 4 arrays as arguments in my odeint and solve the system. For example, at every time step the odeint takes (one from the 8000) to solve the system, I want it to use a different value from these arrays. Is there any way to do it automatically without getting lost in for loops? I tried to do it like this (see code below) but I get the error:

if g2>0: ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

g2 supposed to be 1x1 size at every loop of odeint. So it has to be something with the way I use the 4 arrays (xdot,ydot,xdotdot,ydotdot).

I am new to python and I use python 2.7.12 on Ubuntu 16.04 LTS. Thank you in advance.

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

added_mass_x = 0.03 # kg
added_mass_y = 0.04
mb = 0.3 # kg
m1 = mb-added_mass_x
m2 = mb-added_mass_y
l1 = 0.07 # m
l2 = 0.05 # m
J = 0.00050797 # kgm^2
Sa = 0.0110 # m^2
Cd = 2.44
Cl = 3.41
Kd = 0.000655 # kgm^2
r = 1000 # kg/m^3

c1 = 0.5*r*Sa*Cd
c2 = 0.5*r*Sa*Cl
c3 = 0.5*mb*(l1**2)
c4 = Kd/J
c5 = (1/(2*J))*(l1**2)*mb*l2
c6 = (1/(3*J))*(l1**3)*mb

theta_0 = 10*(np.pi/180) # rad
theta_A = 20*(np.pi/180) # rad
f = 2 # Hz

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

t = np.linspace(0,100,8000) # s

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

# Save data from txt files into numpy arrays
xdot_list = []
ydot_list = []
xdotdot_list = []
ydotdot_list = []

with open('xdot.txt', 'r') as filehandle:
    filecontents = filehandle.readlines()

    for line in filecontents:
        current_place = line[:-1]
        xdot_list.append(current_place)

xdot = np.array(xdot_list, dtype=np.float32) 

with open('ydot.txt', 'r') as filehandle:
    filecontents = filehandle.readlines()

    for line in filecontents:
        current_place = line[:-1]
        ydot_list.append(current_place)

ydot = np.array(ydot_list, dtype=np.float32)

with open('xdotdot.txt', 'r') as filehandle:
    filecontents = filehandle.readlines()

    for line in filecontents:
        current_place = line[:-1]
        xdotdot_list.append(current_place)

xdotdot = np.array(xdotdot_list, dtype=np.float32) 

with open('ydotdot.txt', 'r') as filehandle:
    filecontents = filehandle.readlines()

    for line in filecontents:
        current_place = line[:-1]
        ydotdot_list.append(current_place)

ydotdot = np.array(ydotdot_list, dtype=np.float32) 

def inverse(k,t,xdot,ydot,xdotdot,ydotdot):
    vcx_i = k[0]
    vcy_i = k[1]
    psi_i = k[2]
    wz_i = k[3]
    theta_i = k[4]
    theta_deg_i = k[5]

    # Subsystem 4
    vcx_i = xdot*np.cos(psi_i)-ydot*np.sin(psi_i)
    vcy_i = ydot*np.cos(psi_i)+xdot*np.sin(psi_i)
    psidot_i = wz_i
    vcxdot_i = xdotdot*np.cos(psi_i)-xdot*np.sin(psi_i)*psidot_i-ydotdot*np.sin(psi_i)-ydot*np.cos(psi_i)*psidot_i
    vcydot_i = ydotdot*np.cos(psi_i)-ydot*np.sin(psi_i)*psidot_i+xdotdot*np.sin(psi_i)+xdot*np.cos(psi_i)*psidot_i

    g1 = -(m1/c3)*vcxdot_i+(m2/c3)*vcy_i*wz_i-(c1/c3)*vcx_i*np.sqrt((vcx_i**2)+(vcy_i**2))+(c2/c3)*vcy_i*np.sqrt((vcx_i**2)+(vcy_i**2))*np.arctan2(vcy_i,vcx_i)
    g2 = (m2/c3)*vcydot_i+(m1/c3)*vcx_i*wz_i+(c1/c3)*vcy_i*np.sqrt((vcx_i**2)+(vcy_i**2))+(c2/c3)*vcx_i*np.sqrt((vcx_i**2)+(vcy_i**2))*np.arctan2(vcy_i,vcx_i)

    A = 12*np.sin(2*np.pi*f*t+np.pi) # eksiswsi tail_frequency apo simulink
    if A>=0.1:
        wzdot_i = ((m1-m2)/J)*vcx_i*vcy_i-c4*wz_i**2*np.sign(wz_i)-c5*g2-c6*np.sqrt((g1**2)+(g2**2))
    elif A<-0.1:
        wzdot_i = ((m1-m2)/J)*vcx_i*vcy_i-c4*wz_i**2*np.sign(wz_i)-c5*g2+c6*np.sqrt((g1**2)+(g2**2))
    else:
        wzdot_i = ((m1-m2)/J)*vcx_i*vcy_i-c4*wz_i**2*np.sign(wz_i)-c5*g2

    # Subsystem 5
    if g2>0:
        theta_i = np.arctan2(g1,g2)
    elif g2<0 and g1>=0:
        theta_i = np.arctan2(g1,g2)-np.pi
    elif g2<0 and g1<0:
        theta_i = np.arctan2(g1,g2)+np.pi
    elif g2==0 and g1>0:
        theta_i = -np.pi/2
    elif g2==0 and g1<0:
        theta_i = np.pi/2
    elif g1==0 and g2==0:
        theta_i = 0

    theta_deg_i = (theta_i*180)/np.pi


    return [vcxdot_i, vcydot_i, psidot_i, wzdot_i, theta_i, theta_deg_i]

# initial conditions
vcx_i_0 = 0.1257
vcy_i_0 = 0
psi_i_0 = 0
wz_i_0 = 0
theta_i_0 = 0
theta_deg_i_0 = 0
#theta_i_0 = 0.1745
#theta_deg_i_0 = 9.866
k0 = [vcx_i_0, vcy_i_0, psi_i_0, wz_i_0, theta_i_0, theta_deg_i_0]

# epilysi systimatos diaforikwn
k = odeint(inverse, k0, t, args=(xdot,ydot,xdotdot,ydotdot), tfirst=False)

# apothikeysi apotelesmatwn
vcx_i = k[:,0]
vcy_i = k[:,1]
psi_i = k[:,2]
wz_i = k[:,3]
theta_i = k[:,4]
theta_deg_i = k[:,5]

# Epanalipsi tu Subsystem 5 gia na mporun na plotaristun ta theta_i, theta_deg_i
theta_i = [inverse(k_i, t_i)[4] for t_i, k_i in zip(t, k)]
theta_deg_i = [inverse(k_i, t_i)[5] for t_i, k_i in zip(t, k)]

# Ypologismos mesis gwnias theta kai platus talantwsis
mesi_gwnia = sum(theta_i)/len(theta_i) # rad
platos = (max(theta_i)-min(theta_i))/2

UPDATE: The most relevant solution I found so far is this: Solving a system of odes (with changing constant!) using scipy.integrate.odeint? But since I have only values of my variables in arrays and not the equation of the variables that depend on time (e.g. xdot=f(t)), I tried to aply an interpolation between the values in my arrays, as shown here: ODEINT with multiple parameters (time-dependent) I managed to make the code running without errors, but the total time increased dramatically and the results of the system solved are completely wrong. I tried any possible type of interpolation that I found here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html but still wring outcome. That means that my interpolation isn't the best possible, or my points in the arrays (8000 values) are too much to interpolate between them and solve the system correctly.