Question: I am modelling a Coupled Spring-Mass-System: Two objects with masses m1 and m2, and are coupled through springs with spring constants k1 and k2, with damping d1 and d2. Method 1: taking a cookbook-script using ODEINT to solve the differential equations. # Use ODEINT to solve the differential equations defined by the vector field from scipy.integrate import odeint import matplotlib.pyplot as plt import pandas as pd import numpy as np def vectorfield(w, t, p): """ Defines the differential equations for the coupled spring-mass system. Arguments: w : vector of the state variables: w = [x1,y1,x2,y2] t : time p : vector of the parameters: p = [m1,m2,k1,k2,b1,b2] """ x1, y1, x2, y2 = w m1, m2, k1, k2, b1, b2 = p f = [y1, (-d1 * y1 - k1 * x1 + k2 * (x2 - x1 )+ d2 * (y2-y1)) / m1, y2, (-d2 * (y2-y1) - k2 * (x2 - x1 )) / m2] return f # 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) # Initial conditions # x1 and x2 are the initial displacements; y1 and y2 are the initial velocities x1 = 0.5 y1 = 0.0 x2 = 0.25 y2 = 0.0 # ODE solver parameters abserr = 1.0e-8 relerr = 1.0e-6 stoptime = 100.0 numpoints = 5001 # Create the time samples for the output of the ODE solver. # I use a large number of points, only because I want to make # a plot of the solution that looks nice. t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)] # Pack up the parameters and initial conditions: p = [m1, m2, k1, k2, d1, d2] w0 = [x1, y1, x2, y2] # Call the ODE solver. wsol = odeint(vectorfield, w0, t, args=(p,), atol=abserr, rtol=relerr) # convert zip object to tuple object temp=tuple(zip(t, wsol[:,0],wsol[:,1],wsol[:,2],wsol[:,3])) # convert tulpe object to list then to dataframe df=pd.DataFrame(list(map(list,temp))) # Create plots with pre-defined labels. fig, ax = plt.subplots() ax.plot(df.loc[:,0], df.loc[:,1],label='displacement object1') 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] or V [m/s]', fontdict=None, labelpad=None, loc='center')  Method 2: Using Runge-Kutta-4th Order (programmed by myself) to solve the differential equations 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_sys1(f, x1, y1, x2, y2, h): k1 = f(x1,y1,x2,y2) k2 = f(x1,y1+k1*h/2,x2,y2) k3 = f(x1,y1+k2*h/2,x2,y2) k4 = f(x1,y1+k3*h,x2,y2) temp1= y1 temp2= y1+k1*h/2 temp3= y1+k2*h/2 temp4= y1+k3*h displace_solution = x1 + (temp1 + 2 * temp2 + 2 * temp3 + temp4)*h / 6 velocity_solution = y1 + (k1 + 2 * k2 + 2 * k3 + k4)*h / 6 acc_solution = f(x1,y1,x2,y2) return displace_solution, velocity_solution, acc_solution def runge_kutta_4_sys2(f, x1, y1, x2, y2, h): k1 = f(x1,y1,x2,y2) k2 = f(x1,y1,x2,y2+k1*h/2) k3 = f(x1,y1,x2,y2+k2*h/2) k4 = f(x1,y1,x2,y2+k3*h) temp1= y2 temp2= y2+k1*h/2 temp3= y2+k2*h/2 temp4= y2+k3*h displace_solution = x2 + (temp1 + 2 * temp2 + 2 * temp3 + temp4)*h / 6 velocity_solution = y2 + (k1 + 2 * k2 + 2 * k3 + k4)*h / 6 acc_solution = f(x1,y1,x2,y2) return displace_solution, velocity_solution, acc_solution x1 = 0.5 y1 = 0.0 x2 = 0.25 y2 = 0.0 h = 0.02 numpoints = 5000 time = 0 temp=0 df = pd.DataFrame(index=range(1+numpoints),columns=range(7)) df.iloc[0]=[time,x1,y1,temp,x2,y2,temp] for i in range(numpoints): x1, y1, acc1 = runge_kutta_4_sys1(system1, x1, y1, x2, y2, h) x2, y2, acc2 = runge_kutta_4_sys2(system2, x1, y1, x2, y2, h) time=time+h df.iloc[i+1]=[time,x1,y1,acc1,x2,y2,acc2] df.iloc[i,3]=acc1 df.iloc[i,6]=acc2 # Create plots with pre-defined labels. fig, ax = plt.subplots() ax.plot(df.loc[:,0], df.loc[:,1],label='displacement object1') 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] or V [m/s]', fontdict=None, labelpad=None, loc='center')  There is clearly deviation in the solutions between method 1 and method 2.  I assume the deviation is due to inaccurate Runge-Kutta-integration. Is my assumption correct? How to set or improve the accuracy for the Runge-Kutta-Methode? I tried searching in internet to found the possible restricts and accuracy using runger-kutta method, without success.

Answer: As @lastchance first observed in the comments, your RK4 is not correct as it doesn't apply the method to all the four equations at the same time. I also don't understand the rationale behind the two velocity_solution equations. Here's RK4 method applied to the system of four equations. For pedagogical reasons I add it explicitly, in all its "glorious" complexity. For slightly more complex systems, one would have to add abstractions that would allow the groups of k factors to be computed in loops. 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 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] df.iloc[i, 3] = acc1 df.iloc[i, 6] = acc2