Radio active decay - Monte carlo method

from pylab import*
N0=1000
la=2
dt=0.05
tmax=50
N=N0
t=0
while t<=tmax:
    plot(t,N,'.b')
    for i in range(1,N):
        x=random()
       
        if x<la*dt:
            N=N-1
   
    t=t+dt
td=linspace(0,50,1000)
nd=N0*exp(-la*td)
plot(td, nd, 'r-', linewidth=2)
show()