Periodic boundary condition and 4th order Runge Kutta method used for solving partial differential equation(1D wave packet propagation)

Question

# RK4         
# disp = u0 # initial condition array
flag = True
#h = 0
timeee = 0.000
timemax = 5001 #ms
while(timeee<timemax):
    predisp = np.zeros(N+4)
    ddisp = np.zeros(N+2) # derivative of displacement
    k1 = np.zeros(N)
    K2 = np.zeros(N)
    K3 = np.zeros(N)
    K4 = np.zeros(N)
    # stage 1
    for l in range (0,N):
        predisp[l+2] = u0[l]

    predisp[N+1] = predisp[3]
    predisp[N+2] = predisp[3]
    predisp[N+3] = predisp[4]
    predisp[3] = predisp[N]
    predisp[0] = predisp[N-1]

    for p in range (0,N):
        ddisp[p+1] = uj_prime[p]

    ddisp[N+1] = ddisp[3]
    ddisp[N] = ddisp[3]
    ddisp[0] = ddisp[N+1]     

    for i in range(1,N+1):
        k1[i-1] = -c*ddisp[i]
    k1[N-1] = k1[0]

    for y in range (0,N-1):
        u0[y] = predisp[y+2]+ 0.5*(dt/a_sec)*k1[y]
    u0[N-1] = u0[0]

    #Stage 2
    u_initial_colmat = np.transpose(u0)    
    uj_prime = C @ u_initial_colmat

    for k in range(0,N):
        ddisp[k+1] = uj_prime[k]
    ddisp[N+1] = ddisp[3]
    ddisp[N] = ddisp[3]
    ddisp[0] = ddisp[N+1]    

    for i in range(1,N+1):
        K2[i-1] = -c*ddisp[i]
    K2[N-1] = K2[0]
    

    for y in range(0,N-1):
        u0[y] = predisp[y+2] + 0.5*(dt/a_sec)*K2[y]
    u0[N-1] = u0[0]    

    #Stage 3
    u_initial_colmat = np.transpose(u0)    
    uj_prime = C @ u_initial_colmat

    for k in range(0,N):
        ddisp[k+1] = uj_prime[k]
    ddisp[N+1] = ddisp[3]
    ddisp[N] = ddisp[3]
    ddisp[0] = ddisp[N+1] 

    for i in range(1,N+1):
        K3[i-1] = -c*ddisp[i]
    K3[N-1] = K3[0]   


    for y in range(0,N-1):
        u0[y] = predisp[y+2] + 1.0*(dt/a_sec)*K3[y]
    u0[N-1] = u0[0]
    #Stage 4
    u_initial_colmat = np.transpose(u0)    
    uj_prime = C @ u_initial_colmat

    for k in range(0,N):
        ddisp[k+1] = uj_prime[k]
    ddisp[N+1] = ddisp[3]
    ddisp[N] = ddisp[3]
    ddisp[0] = ddisp[N+1]    

    for i in range(1,N+1):
        K4[i-1] = -c*ddisp[i]
    K4[N-1] = K4[0]
    
    for y in range(0,N-1):
        u0[y] = predisp[y+2] + (1/6)*(dt/a_sec)*(k1[y] + 2*K2[y] + 2*K3[y] + K4[y])  
    u0[N-1] = u0[0]

    milli_time = timeee/a_sec        
    uu =  ((np.exp(alpha*((x-c*milli_time-x0)**2))*(np.cos(k0*(x-c*milli_time-x0))))) # Analytical Solution
    timeee = timeee + dt
    if (timeee >= timemax):
            exit()

It is the part of the code where RK4 is used. Tangent hyperbolic function(double-sided) is used as the grid. I need to solve the 1D wave equation using Runge Kutta of order 4 for time discretisation with periodic boundary conditions. I know that while getting matrices for solving derivatives is correct and the periodic boundary conditions are applied properly. But during the application of the RK4 method, I am a bit unsure if the periodic conditions are applied properly. spatial discretisation algo have derivative terms considering three points on RHS: j-1, j and j+1
N is the number of nodes in the grid. The above image shows how the RK4 method is used in the code.
The derivative term is written as ddisp in the code. Any help is highly appreciated.