plt.gcf().clear()

N = 60
X = np.linspace(-4, 7, N)
Y = np.linspace(-4, 6, N)
X, Y = np.meshgrid(X, Y)

pos = np.empty(X.shape + (2,))
pos[:, :, 0] = X
pos[:, :, 1] = Y

def multivariate_gaussian(pos, mu, Sigma):
    n = mu.shape[0]
    Sigma_det = np.linalg.det(Sigma)
    Sigma_inv = np.linalg.inv(Sigma)
    N = np.sqrt((2*np.pi)**n * Sigma_det)
    fac = np.einsum('...k,kl,...l->...', pos-mu, Sigma_inv, pos-mu)
    return np.exp(-fac / 2) / N

Z1 = multivariate_gaussian(pos, mean1, sd1)
Z2 = multivariate_gaussian(pos, mean2, sd2)
Z = Z2 - Z1

plt.contourf(X, Y, Z)
ax = plt.axes()
ax.arrow(-7*max_w_x, -7*max_w_y, 15*max_w_x, 15*max_w_y, fc='k', ec='k')
plt.show()