import matplotlib.pyplot as plt
   import numpy as np
   from scipy.integrate import odeint
   # Numerical solution using an ordinariy differential equation solver
   def ode_system(initial_cs, t, k1, k2):
       c1, c2, c3, c4, c5 = initial_cs
       dc1_dt = -2*k1*c1**2*c2
       dc2_dt = -k1*c1**2*c2
       dc3_dt = k1*c1**2*c2-2*k2*c3**2
       dc4_dt = k2*c3**2
       dc5_dt = k2*c3**2
       return dc1_dt, dc2_dt, dc3_dt, dc4_dt, dc5_dt
   initial_cs = (0.2, 0.3, 0.0001, 0.0001, 0.0001)
   k1 = 0.01
   k2 = 1e-3
   ts = np.linspace(0, 6000, num=200)
   solution = odeint(ode_system, initial_cs, ts, args=(k1, k2))
   concentration_names = {1:  r"$H_2$", 2: r"$O_2$", 3: r"$H_2O$", 4: r"$H_3O^+$", 5: r"$OH^-$"}
   for i in range(1, 6):
       plt.plot(ts, solution[:, i-1], label=concentration_names[i])
   plt.xlabel(r"time [$\tau$]")
   plt.ylabel("cocentrations")
   plt.title(r"Solution of differential equation system, $k_1=0.01 \frac{\tau^-1}{c_0^2}, k_2=10^{-3} \frac{\tau^-1}{c_0}$")
   plt.legend()
   plt.show()