153556

import numpy as np
import matplotlib.pyplot as plt

Q = 200
P = np.zeros(Q)
P[0] = 1
PP = np.zeros(Q)
t_max = 100
dt = 0.001
lmbd = 4
mu = 5
t = 0
p0 = [1]
p1 = [0]

while t < t_max:
    PP[0] = dt * (-lmbd * P[0] + mu * P[1])
    for j in range(1, Q-1):
        PP[j] = dt * (-(lmbd + mu) * P[j] + lmbd * P[j-1] + mu * P[j + 1])
    for j in range(Q):
        P[j] += PP[j]
        t += dt
    p0.append(P[0])
    p1.append(P[1])
print(P, p0, p1, sep='\n')
t_to_plot = np.linspace(0, 100, 200)
p0 = p0[:200]
p1 = p1[:200]
#plt.plot(t_to_plot, P)
plt.plot(t_to_plot, p0)
plt.plot(t_to_plot, p1)
plt.show()