This is a homework of mine to plot the gaussian and poisson distribution for r = 20, mu = 20 and sigma = sqrt(20). But the problem is that the value when computing of the computer is capped somewhere, and the large number like 20^20 cannot be computed, so the graph looks horribly wrong after r >= 8. Here is my code:

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from scipy.special import factorial

#c)
sigma = np.sqrt(20)
mu1 = 20
b = np.sqrt(2*np.pi)
def fGauss(r):
    return 1/(b * sigma) * np.exp(-(r - mu1)**2 / (2*sigma**2))
def fPoisson2(r):
    return mu1**r * np.exp(-mu1)/factorial(r)

x1 = np.arange(1, 21, 1)
Gauss = [fGauss(r) for r in x1]
Poisson2 = [fPoisson2(r) for r in x1]


plt.bar(x1, Poisson2, color = 'b', alpha = 1, label = 'Poisson-Verteilung')
plt.bar(x1, Gauss, color = 'r', alpha = 0.2, label = 'Gauss-Verteilung')
plt.xlim(-0.1, 20.1)
plt.xticks(x1)
plt.xlabel('r')
plt.ylabel('Wahrscheinlichkeit')
plt.title(r'Poisson- und Gauss-Verteilung für $(mu = 20; r = 1,2,...,20)$')
plt.show()

I have asked a friend of mine for advice and he said to use np.random.poisson, but I still don’t know how to work around it, so I post this here, hoping for help. In attachment is my file.
If anyone has question regarding the task, I will provide as much as I could.
Many thanks