The splitting-operator method has been widely used in the numerical simulation of complex combustion system with its high computation efficiency. The main difficult encountered in the method is how to integrate the stiff chemical kinetics ordinary differential equation (ODE) efficiently and accurately in each grid node of flow field. Although several ODE software packages such as LSODE and VODE have very predominant performance in integrating stiff and non-stiff case, it couldn’t ensure the positivity of mass fraction in integrating stiff chemical kinetic ODE. A detail liquid rocket hypergolic bipropellant combination of Monomethylhydrazine (MMH)/Red Fuming Nitric Acid (RFNA) chemical mechanism consisting of 550 elementary reactions among 83 species was constituted basing on several relevant literatures firstly. And then a modified implicit iterative difference algorithm with variable time steps was developed to tackle the problem of mass fraction non-positive. The results from the integration of MMH /nitrogen tetroxide (NTO) chemical kinetic ODE indicate that the algorithm is always retain stability and positive value of mass fraction. The time step length would also be automatically adjusted at different chemical reaction time in the algorithm to keep its computational efficiency.