1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
# Radio active decay. Solve 2 equations DN/dt = -L * N at a time
from pylab import *
from scipy import integrate

L1 = .5          # decay constant
N1 = 1000        # value at t = 0
L2 = 1.0
N2 = 2000

def derivative(Y, t0):
    return [-L1 * Y[0], -L2 * Y[1] ]  # dN/dt = -L * N , radioactive decay

t = arange(0, 3, 0.1)                      # time span and steps
nt0 = [N1, N2]                             # values of Ns at t=0
nt = integrate.odeint(derivative, nt0, t)  # integrate

plot(t, nt[:,0])  # extract the first column from the 2D array
plot(t, nt[:,1])  # extract the second column from the 2D array
show()