System of differential equations with time depende

Let's say that I have a system of differential equations and I want to solve it with odeint. Some of the system's constants are time depended and I have their values stored in arrays (a,b,c and d with shape (8000, ) ). I want the system to use a different value of these constants at each time step. See the simplified code example:

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

a = np.array([0, 5, 6, 12, 1.254, ..., 0.145])     # shape (8000, )
b = np.array([1.45, 5.9125, 1.367, ..., 3.1458])
c = np.array([0.124, 0.258, 0.369, ..., 0.147])
d = np.array([7.145, 5.123, 6.321, ..., 0.125])

def system(k,t):
    vcx_i = k[0]
    vcy_i = k[1]
    psi_i = k[2]
    wz_i = k[3]

    vcx_i_dot = a*np.cos(psi_i)-b*np.sin(psi_i)
    vcy_i_dot = b*np.cos(psi_i)+a*np.sin(psi_i)
    psi_i_dot = wz_i
    wz_i_dot = c*vcx_i-a*np.sin(psi_i)-d*np.sin(psi_i)-b*np.cos(psi_i)

    return [vcx_i_dot, vcy_i_dot, psi_i_dot wz_i_dot]

k0 = [0.1257, 0, 0, 0]

k = odeint(system, k0, t)

vcx_i = k[:,0]
vcy_i = k[:,1]
psi_i = k[:,2]
wz_i = k[:,3]

psi_i = [system(t_i, k_i)[2] for k_i, t_i in zip(k,t)]
wz_i = [system(t_i, k_i)[3] for k_i, t_i in zip(k,t)]

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. a=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 wrong 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. What is the right way to solve a system like this? I am new to python and I use python 2.7.12 on Ubuntu 16.04 LTS. Thank you in advance.

The interpolators are usually very fast, so probably there is something else in your function. You can however, attempt different interpolators (like InterpolatedUnivariateSpline), or reducing the interpolation nodes to increase speed. But I would aim at your integration instead.

Lately, ode and odeint have been replaced by other, more flexible functions (see here)

I would start with and explicit method instead of an implicit one (default for solve_ivp is runge kutta, while default for odeint is LSODA):

interp = scipy.interpolate.interp1d(t,(a,b,c,d))

def system(t,k):
    vcx_i = k[0]
    vcy_i = k[1]
    psi_i = k[2]
    wz_i = k[3]
    a, b, c, d = interp(t)

    vcx_i_dot = a*np.cos(psi_i)-b*np.sin(psi_i)
    vcy_i_dot = b*np.cos(psi_i)+a*np.sin(psi_i)
    psi_i_dot = wz_i
    wz_i_dot = c*vcx_i-a*np.sin(psi_i)-d*np.sin(psi_i)-b*np.cos(psi_i)
    return [vcx_i_dot, vcy_i_dot, psi_i_dot, wz_i_dot]

k0 = [0.1257, 0, 0, 0]
steps = 1
method = 'RK23'
atol = 1e-3
s = solve_ivp(dydt, (0, 100), k0, method=method, t_eval=t, atol=atol, vectorized=True)

you can increase / reduce atol or change the method.