2D lattice with fixed boundary condition
Testing script from a few days ago:

import numpy as  np


original:

def CML2d(x):
    eps = 0.3
    r = 3.9
    xn = np.zeros(N+1, float)
    for i in range(1, N+1):
        if i&gt;L and i&lt;=(L-1)*L:
            if i%L==1:
                xl, xr = 0., x[i+1]
                xu, xd = x[i-L], x[i+L]
            elif i%L==0:
                xl, xr = x[i-1], 0.
                xu, xd = x[i-L], x[i+L]
            else:
                xl, xr = x[i-1], x[i+1]
                xu, xd = x[i-L], x[i+L]
        elif i&gt;1 and i&lt;L:
            xl, xr = x[i-1], x[i+1]
            xu, xd = 0., x[i+L]
        elif i&gt;(L-1)*L+1 and i&lt;L*L:
            xl, xr = x[i-1], x[i+1]
            xu, xd = x[i-L], 0.
        elif i==1:
            xl, xr = 0., x[i+1]
            xu, xd = 0., x[i+L]
        elif i==L:
            xl, xr = x[i-1], 0.
            xu, xd = 0., x[i+L]
        elif i==(L-1)*L+1:
            xl, xr = 0., x[i+1]
            xu, xd = x[i-L], 0.
        elif i==L*L:
            xl, xr = x[i-1], 0.
            xu, xd = x[i-L], 0.
        value = (1-eps)*r*x[i]*(1-x[i]) + 0.25*eps*( r*xl*(1-xl) + r*xr*(1-xr) + r*xu*(1-xu) + r*xd*(1-xd) )
        xn[i] = value
    return xn


partial attempt to work with 2d x; the intention was to replace all the uses of L with 2d array indexing tests.

def CML2d_1(x):
    eps = 0.3
    r = 3.9
    L,_ = x.shape
    N = L*L
    # xn = np.zeros((L,L), float)
    xn = (1-eps)*r*x*(1-x)
    x = x.flat
    #for i in range(N):
    for j in range(L):
      for k in range(L):
        i = k+L*j
        i1 = i+1
        if i1&gt;L and i1&lt;=(L-1)*L:
            if i1%L==1:
                xl, xr = 0., x[i+1]
                xu, xd = x[i-L], x[i+L]
        ....
        #value = (1-eps)*r*x[i]*(1-x[i])
        value = 0.25*eps*( r*xl*(1-xl) + 
                            r*xr*(1-xr) + 
                            r*xu*(1-xu) + 
                            r*xd*(1-xd) )
        #xn.flat[i] = value
        xn[j,k] += value
    return xn


But I then realized that I don't need to iterate over all points.  Instead I could just sum subarrays  (similar to the array of taking a 1d diff, x[1:]-x[:-1]:

# x[i+1] = (1-eps)*   r*x[i]*(1-x[i]) + 
#          0.25*eps*( r*x[i-1]*(1-x[i-1]) + 
#                     r*x[i+1]*(1-x[i+1]) + 
#                     r*x[i-L]*(1-x[i-L]) + 
#                     r*x[i+L]*(1-x[i+L]) )

def CML2d_2(x):
    eps = 0.3
    r = 3.9
    x2 = r * x * (1-x)
    xn = (1-eps) * x2
    xn[1:,:]  += 0.25 * eps * x2[:-1,:]
    xn[:-1,:] += 0.25 * eps * x2[1:,:]
    xn[:,1:]  += 0.25 * eps * x2[:,:-1]
    xn[:,:-1] += 0.25 * eps * x2[:,1:]
    return xn


then making use of a covolve2d function in scipy.signal:

from scipy import signal
def CML2d_3(x):
    eps = 0.3
    r = 3.9
    in2 = np.zeros((3,3))
    in2[1,:] = 0.25 * eps
    in2[:,1] = 0.25 * eps
    in2[1,1] = (1-eps)
    print(in2)
    x2 = r * x * (1-x)
    res = signal.convolve2d(x2, in2, mode='same', boundary='fill', fillvalue=0)
    return res  


testing:

L = 10 #side of 2d lattice
#L = 4
N = L*L #number of sites in 2d lattice
x0 = np.random.uniform(0.1, 0.9, N+1) #initial values for x
x0[0]=np.nan  # test if this is used
res=CML2d(x0.tolist())
print(res)

x2d = x0[1:].reshape(L,L)
res1=CML2d_1(x2d)
print(x0.shape)
print(res1.shape)
print(np.allclose(res1.flatten(), res[1:]))

print(np.allclose(res1, CML2d_2(x2d)))
print(np.allclose(res1, CML2d_3(x2d)))