zn+1 = zn*r*(1-zn)

population.py

%matplotlib inline

import matplotlib
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib

L = 5000
N = 70
Ζ0 = 0.0000001
fig=plt.figure(figsize=(18, 16), dpi= 80, facecolor='w', edgecolor='k')
Rs = np.linspace(1, 4., L)
Zn = np.ones(L) * Ζ0
rainbow = cm.rainbow(np.linspace(0, 1, N))
for c in rainbow:
    Zn = Rs * Zn * (1 - Zn)
    ax = plt.scatter(Rs, Zn , s=0.05, color=c)
cb = plt.colorbar(cmap=rainbow, norm=matplotlib.colors.Normalize(vmin=0, vmax=1/N), orientation='horizontal')
cb.set_label("Iteration # ")