were firstly solved using a method that automatically switches between nonstiff (Adams) and stiff (BDF) solvers, as implemented in scipy.integrate.odeint in Python. The steady-state coverages were further converged with a system of algebraic equations using the coverages solved from ODE as initial guesses. Newton's method implemented in mpmath.findroot in Python was used.