Multivariate Normal Distribution Function Matrix multiplication behaving not as expected

So I’ve create a function that aims to generate a 2D Gaussian array, as per the equation,
$f_X(x_1, x_2,…, x_k)=\frac{exp(-\frac{1}{2}}(x-\mu)^T\Sigma^{-1}(x-\mu)){\sqrt{(2\pi)^k|\Sigma|}}$

However, the peaks are split at the edges of the parameter space instead of where the mean is.

def Gauss_2D(x_grid, y_grid, mean=np.array((0, 0)), covariance_matrix=np.array([[1, 0], [0, 1]])):
    x=np.linspace(-3, 3, x_grid)-mean[0]
    y=np.linspace(-3, 3, y_grid)-mean[1]

    xy=np.vstack((x, y))
    inv_cov_matrix = np.linalg.inv(covariance_matrix)

    exponent = -0.5 * (xy.T @ inv_cov_matrix @ xy)
    z = (1 / (2 * np.pi * np.sqrt(np.linalg.det(covariance_matrix)))) * np.exp(exponent)

    return z

I would also like to mention that I do have a code that works, shown below for posterity sake, I just cannot seem to understand why the code above doesn’t work.

def Gauss2D(x_grid, y_grid, mean=(0, 0), covariance_matrix=[[1, 0], [0, 1]]):
    x = np.linspace(-3, 3, x_grid)
    y = np.linspace(-3, 3, y_grid)
    x, y = np.meshgrid(x, y)

    xy = np.column_stack((x.flatten(), y.flatten()))
    inv_cov_matrix = np.linalg.inv(covariance_matrix)
    
    exponent = -0.5 * np.sum((xy - mean) @ inv_cov_matrix * (xy - mean), axis=1)
    z = (1 / (2 * np.pi * np.sqrt(np.linalg.det(covariance_matrix)))) * np.exp(exponent)

    return z.reshape(rows, cols)

The output of,

g2d= Gauss_2D(100, 100)
plt.imshow(g2d, cmap='viridis', extent=(-3, 3, -3, 3))

yields a gaussian centered on the bottom left and top right corners instead of the mean location when using the 1st function.