What is wrong with this Fourier transform? (in python)

As people who have read my previous posts on this site could tell, I have been trying to implement a PDE solver that uses FFT in Python. The programming part is mostly worked out, but the program produces an (very suitable for this site) overflow error (basically it grows very much until it becomes a NaN).

After ruling out all other possibilities, I pinned down the problem to the FFT and the way I am trying to do the derivatives, so I decided to test two different FFT's (numpy's fft module and the pyFFTW package) with the following code:

import pyfftw
import numpy as np
import matplotlib.pyplot as plt


def fftw_(y: np.ndarray) -> np.ndarray:
    a = pyfftw.empty_aligned((N, N), dtype=np.float64)
    b = pyfftw.empty_aligned((N, N//2+1), dtype=np.complex128)
    fft_object = pyfftw.FFTW(a, b, axes=(0, 1), direction='FFTW_FORWARD', flags=('FFTW_MEASURE',), threads=12)
    y_hat = fft_object(y)
    return y_hat


def ifftw_(y_hat: np.ndarray) -> np.ndarray:
    a = pyfftw.empty_aligned((N, N//2+1), dtype=np.complex128)
    b = pyfftw.empty_aligned((N, N), dtype=np.float64)
    fft_object = pyfftw.FFTW(a, b, axes=(0, 1), direction='FFTW_BACKWARD', flags=('FFTW_MEASURE',), threads=12)
    y = fft_object(y_hat)
    return y


def func(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    return np.exp(x)*np.sin(y)

dx = 0.02

x = np.arange(-1, 1, dx)
y = np.arange(-1, 1, dx)
X, Y = np.meshgrid(x, y)

N = len(x)

kxw, kyw = np.meshgrid(np.fft.rfftfreq(N, dx), np.fft.fftfreq(N, dx))
Lapw = -4*np.pi**2*(kxw**2+kyw**2)

kxnp, kynp = np.meshgrid(np.fft.fftfreq(N, dx), np.fft.fftfreq(N, dx))
Lapnp = -4*np.pi**2*(kxnp**2+kynp**2)



z = func(X, Y)

lap_z_w = ifftw_(Lapw*fftw_(z))
lap_z_np = np.fft.ifft2(Lapnp*np.fft.fft2(z))

lap_z_np_mag = np.abs(lap_z_np)
lap_z_np_ang = np.angle(lap_z_np)

plt.imshow(z, cmap='plasma')
plt.colorbar()
plt.savefig("f.png", dpi=200)
plt.clf()

plt.imshow(lap_z_w, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_fftw.png", dpi=200)
plt.clf()

plt.imshow(lap_z_np_mag, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_np_mag.png", dpi=200)
plt.clf()

plt.imshow(lap_z_np_ang, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_np_ang.png", dpi=200)
plt.clf()

Here the np.ndarray's named Lapw and Lapnp are what I thought should do the discrete Laplacian. And the function I chose, eˣsin(y), is a harmonic function, so its Laplacian should be zero.

But the results from the program are very far from this expected value. In Particular I get:

The original function f

The "Laplacian" of f with pyFFTW

The magnitude and phase of the "Laplacian" of f with Numpy

Looking at the values of these plots (please do note the range in the colorbar and the fact that 20000 is not any kind of decent approximation to 0) makes it clear why the program I made is giving an overflow, but I don't know how to correct this. Any help would be greatly appreciated.

The mistake here is that your assumptions about your function are not correct. e^x sin(y) might seem harmonic, but you only calculated it for -1 < x,y < 1. The fft will implicitly continue it periodically, i.e. you get discontinuities at all the edges of your function. If the function is not continuous, it's not harmonic and especially you get divergences in the Fourier transfom. This is what makes your FFT diverge at the edges. Besides that, "far away" from the edges, the results look as expected.