import matplotlib.pyplot as plt
from photutils.datasets import make_noise_image
from photutils.psf import CircularGaussianPRF, make_psf_model_image

shape = (256, 256)
fwhm = 4.7
psf_model = CircularGaussianPRF(fwhm=fwhm)
psf_shape = (11, 11)
n_sources = 100
flux = (500, 1000)
border_size = (7, 7)
data, stars = make_psf_model_image(shape, psf_model, n_sources,
                                   flux=flux,
                                   model_shape=psf_shape,
                                   border_size=border_size, seed=123)
noise = make_noise_image(shape, mean=0, stddev=2, seed=123)
data += noise

fig, ax = plt.subplots(figsize=(8, 8))
ax.imshow(data, origin='lower')