Thanks to "Lastchance" for your hint. And Thanks to "Kikon" for the detailed explanation. @Kikon: I have modified your script to see the difference between synchronous and sequential RK4 integration of system 1 &2.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib.pyplot as plt




# Parameter values
tau=1.02 
f_1=0.16 
f_2=tau*f_1

m1 = 2000000 # mass1
m2 = 20000 # mass2

# Spring constants
k1 = m1*pow(2*np.pi*f_1,2)
k2 = m2*pow((2*np.pi*f_2),2)


# damping
d1 = (0.04/2/np.pi)*2*pow(k1*m1,0.5)



d_d1=6000
l_p=9.81/pow(2*np.pi*f_2,2)
b=0.3
d2=d_d1*pow((l_p-b)/l_p,2)

def system1(x1, y1, x2, y2):
    return (-d1 * y1 - k1 * x1 + k2 * (x2 - x1) + d2 * (y2 - y1)) / m1

def system2(x1, y1, x2, y2):
    return (-d2 * (y2 - y1) - k2 * (x2 - x1)) / m2

def runge_kutta_4(f1, f2, x1, y1, x2, y2, h):
    k1x1 = y1
    k1x2 = y2
    k1y1 = f1(x1, y1, x2, y2)
    k1y2 = f2(x1, y1, x2, y2)

    k2x1 = y1 + h * k1y1 / 2
    k2x2 = y2 + h * k1y2 / 2
    k2y1 = f1(x1 + h * k1x1 / 2, y1 + h * k1y1 / 2, x2 + h * k1x2 / 2, y2 + h * k1y2 / 2)
    k2y2 = f2(x1 + h * k1x1 / 2, y1 + h * k1y1 / 2, x2 + h * k1x2 / 2, y2 + h * k1y2 / 2)

    k3x1 = y1 + h * k2y1 / 2
    k3x2 = y2 + h * k2y2 / 2
    k3y1 = f1(x1 + h * k2x1 / 2, y1 + h * k2y1 / 2, x2 + h * k2x2 / 2, y2 + h * k2y2 / 2)
    k3y2 = f2(x1 + h * k2x1 / 2, y1 + h * k2y1 / 2, x2 + h * k2x2 / 2, y2 + h * k2y2 / 2)

    k4x1 = y1 + h * k3y1
    k4x2 = y2 + h * k3y2
    k4y1 = f1(x1 + h * k3x1, y1 + h * k3y1, x2 + h * k3x2, y2 + h * k3y2)
    k4y2 = f2(x1 + h * k3x1, y1 + h * k3y1, x2 + h * k3x2, y2 + h * k3y2)

    x1_next = x1 + h * (k1x1 + 2 * k2x1 + 2 * k3x1 + k4x1) / 6
    x2_next = x2 + h * (k1x2 + 2 * k2x2 + 2 * k3x2 + k4x2) / 6
    y1_next = y1 + h * (k1y1 + 2 * k2y1 + 2 * k3y1 + k4y1) / 6
    y2_next = y2 + h * (k1y2 + 2 * k2y2 + 2 * k3y2 + k4y2) / 6

    acc1_next = f1(x1_next, y1_next, x2_next, y2_next)
    acc2_next = f2(x1_next, y1_next, x2_next, y2_next)

    return x1_next, x2_next, y1_next, y2_next, acc1_next, acc2_next

# runge kutta 4th integration, 
'''
x1: system 1 displacement
y1: system1 velocity
acc1: system 1 acceleration
x2: system 2 displacement
y2: system2 velocity
acc2: system 2 acceleration

h: time step

'''
x1 = 0.5
y1 = 0.0
x2 = 0.25
y2 = 0.0

h = 0.02
numpoints = 5000
time = 0

temp1 = system1(x1, y1, x2, y2)
temp2 = system1(x1, y1, x2, y2)

df = pd.DataFrame(index=range(1 + numpoints), columns=range(7))

df.iloc[0] = [time, x1, y1, temp1, x2, y2, temp2]

for i in range(numpoints):
    x1, x2, y1, y2, acc1, acc2 = runge_kutta_4(system1, system2, x1, y1, x2, y2, h)
    time = time + h
    df.iloc[i + 1] = [time, x1, y1, acc1, x2, y2, acc2]

# runge kutta 4th integration - test with sequential integration of system 1 and system 2
x1 = 0.5
y1 = 0.0
x2 = 0.25
y2 = 0.0

h = 0.02
numpoints = 5000
time = 0

temp1 = system1(x1, y1, x2, y2)
temp2 = system1(x1, y1, x2, y2)

df_seq = pd.DataFrame(index=range(1 + numpoints), columns=range(7))

df_seq.iloc[0] = [time, x1, y1, temp1, x2, y2, temp2]

for i in range(numpoints):
    x1, temp_x2, y1, temp_y2, acc1, temp_acc2 = runge_kutta_4(system1, system2, x1, y1, x2, y2, h)
    temp_x1, x2, temp_y1, y2, temp_acc1, acc2 = runge_kutta_4(system1, system2, x1, y1, x2, y2, h)
    time = time + h
    df_seq.iloc[i + 1] = [time, x1, y1, acc1, x2, y2, acc2]


# Create plots with pre-defined labels.
fig, ax = plt.subplots()
ax.plot(df.loc[:,0], df.loc[:,1],label='displacement TT_synch') 
ax.plot(df_seq.loc[:,0], df.loc[:,1],'--',label='displacement TT_seq') 


legend = ax.legend(loc='upper right', shadow=None, fontsize='small')
ax.set_xlabel('time [s]', fontdict=None, labelpad=None, loc='center')
ax.set_ylabel('pos [m]', fontdict=None, labelpad=None, loc='center')    

the result shows that there is no difference between synchron and sequential calculation. synchron vs sequential

Therefore I have got below conclusions:

1. the deviation I asked at the beginning ("odeint vs runge-kutta-4th") is due to incorrect differential equations
2. the synchronous and sequential calculation for system 1 &2 using RK4 have no impact on calculation results.

What is your opinion?