def _pollution_rhs(t, y):
k = np.array([
0.35e0, 0.266e2, 0.123e5, 0.86e-3, 0.82e-3,
0.15e5, 0.13e-3, 0.24e5, 0.165e5, 0.9e4,
0.22e-1, 0.12e5, 0.188e1, 0.163e5, 0.48e7,
0.35e-3, 0.175e-1, 0.1e9, 0.444e12, 0.124e4,
0.21e1, 0.578e1, 0.474e-1, 0.178e4, 0.312e1,
])
r = np.zeros(25)
r[0] = k[0] * y[0]
r[1] = k[1] * y[1] * y[3]
r[2] = k[2] * y[4] * y[1]
r[3] = k[3] * y[6]
r[4] = k[4] * y[6]
r[5] = k[5] * y[6] * y[5]
r[6] = k[6] * y[8]
r[7] = k[7] * y[8] * y[5]
r[8] = k[8] * y[10] * y[1]
r[9] = k[9] * y[10] * y[0]
r[10] = k[10] * y[12]
r[11] = k[11] * y[9] * y[1]
r[12] = k[12] * y[13]
r[13] = k[13] * y[0] * y[5]
r[14] = k[14] * y[2]
r[15] = k[15] * y[3]
r[16] = k[16] * y[3]
r[17] = k[17] * y[15]
r[18] = k[18] * y[15]
r[19] = k[19] * y[16] * y[5]
r[20] = k[20] * y[18]
r[21] = k[21] * y[18]
r[22] = k[22] * y[0] * y[3]
r[23] = k[23] * y[18] * y[0]
r[24] = k[24] * y[19]
dy = np.zeros(20)
dy[0] = -r[0] - r[9] - r[13] - r[22] - r[23] + r[1] + r[2] + r[8] + r[10] + r[11] + r[21] + r[24]
dy[1] = -r[1] - r[2] - r[8] - r[11] + r[0] + r[20]
dy[2] = -r[14] + r[0] + r[16] + r[18] + r[21]
dy[3] = -r[1] - r[15] - r[16] - r[22] + r[14]
dy[4] = -r[2] + 2*r[3] + r[5] + r[6] + r[12] + r[19]
dy[5] = -r[5] - r[7] - r[13] - r[19] + r[2] + 2*r[17]
dy[6] = -r[3] - r[4] - r[5] + r[12]
dy[7] = r[3] + r[4] + r[5] + r[6]
dy[8] = -r[6] - r[7]
dy[9] = -r[11] + r[6] + r[8]
dy[10] = -r[8] - r[9] + r[7] + r[10]
dy[11] = r[8]
dy[12] = -r[10] - r[11] + r[9]
dy[13] = -r[12] + r[11]
dy[14] = r[12]
dy[15] = -r[17] - r[18] + r[15]
dy[16] = -r[19] + r[17]
dy[17] = r[19]
dy[18] = -r[20] - r[21] - r[23] + r[22] + r[24]
dy[19] = -r[24] + r[23]
return dy