Abstract: Numerical modeling of nonequilibrium state-to-state carbon dioxide kinetics is a challenging time-consuming computational task that involves solving a huge system of stiff differential equations and requires optimized methods to solve it. In the present study, we propose and analyse optimizations for the Extended Backward Differential Formula (EBDF) scheme. Using adaptive timesteps instead of fixed ones reduces the number of steps in the algorithm many thousands of times, although with an increase in step complexity. The use of parallel computations to calculate relaxation terms allows one to further reduce the computation time. Numerical experiments on the modeling of spatially homogeneous carbon dioxide vibrational relaxation were performed for optimized computational schemes of different orders. Based on them, the most optimal algorithm of calculations was recommended: a parallel EBDF scheme of fourth-order with an adaptive timestep. This method takes less computational time and memory costs and has the high stability.