Расчет лапласиана функции через метод RBF-FD

Я пытаюсь найти лапласиан функции, используя метод RBF-FD на python. Целевая функция sin(pi*x)*cos(pi*z/2), аналитически найденный лапласиан -5/4*(pi^2)*sin(pi*x)*cos(pi*z/2). Исходные данные:

nx = 101         
nz = nx
dx = 0.01
dz = dx
x    = np.linspace(0, 1, nx)
z    = np.linspace(0, 1, nz)
d2px = np.zeros((nz, nx)) 
d2pz = np.zeros((nz, nx))
X, Z = np.meshgrid(x, z)
points = np.column_stack([X.ravel(), Z.ravel()])
stencil = 5
size = len(points)
def func(x,z):
    return np.sin(np.pi*x)*np.cos(np.pi*z/2)
f = func(points[:, 0], points[:, 1]).reshape(nx, nz)
def laplace_abs(x, z):
    return -5/4*(np.pi**2)*np.sin(np.pi*x)*np.cos(np.pi*z/2)
laplace_abs = laplace_abs(points[:, 0], points[:, 1]).reshape(nx, nz)

Затем я пытаюсь найти коэффициенты RBF-FD на 5-точечном гриде, как в FDM. После нахождения коэффициентов значение лапласиана должно быть представлено в виде взвешенной суммы. Соответствующие формулы представлены здесь:

Основная часть кода с вычислением лапласиана через RBF-FD

def get_fd_indices():
    nbrs = NearestNeighbors(n_neighbors=stencil, algorithm='ball_tree', metric='euclidean').fit(points)
    return nbrs.kneighbors(points)
distances, indices = get_fd_indices()
r = distances
RHS = np.where(distances != 0, r**2*(16*np.log(r)+8), 0)
def LHS(x, xi):
    R = np.zeros((size,stencil,stencil))
    for i in range(size):
        R[i] = distance.cdist(x[i], x[i], 'euclidean')
    LHS = np.where(R != 0, R**4*np.log(R), 0)
    return LHS
LHS = LHS(points[indices], points[indices])
def get_coef(LHS, RHS):
    c = np.zeros((size, stencil, 1))
    for i in range(size):
        Left = LHS[i]
        Right = RHS[i].reshape((stencil,1))
        c[i] = LA.pinv(Left).dot(Right)
    return c
c = get_coef(LHS, RHS)
values = f.ravel()
laplaceRBF = np.zeros(size) 
for i in range(size):
    index = indices[i]
    laplaceRBF[i] = np.sum(values[index]*c[i])
laplaceRBF = laplaceRBF.reshape(nx, nz)

Ожидается, что конечный результат должен быть похож на FDM

def FDM(f, dx):
    for i in range(1, nx - 1):
        d2px[i, :] = (f[i - 1, :] - 2 * f[i, :] + f[i + 1, :]) / dx ** 2 
    for j in range(1, nz - 1):
        d2pz[:, j] = (f[:, j - 1] - 2 * f[:, j] + f[:, j + 1]) / dz ** 2 
    return d2px+d2pz
FDM = FDM(f, dx)

На данный момент результаты полученный по RBF-FD и FDM не сходятся. Подскажите, где может быть ошибка.