Here is a working code and the result shown below, just analyze it. Note that the normalization for the probability density function $p(x)$ used here is $\int p(x) dx$=1.
import numpy as np
import matplotlib.pyplot as plt
n, bins, patches = plt.hist(x=vals, bins='auto', density=True, color='#0504aa', alpha=0.7, rwidth=0.85)
What exactly do you think the formula
does, especially the denominator? For the correct physics it should be the third power of the Euclidean distance.