Local thermodynamic equilibrium is a foundational concept in plasma physics and heat transfer, describing a state where each small region of a system can be treated as if it is in thermodynamic equilibrium, even if the whole system is not. However, achieving accurately perfect thermodynamic equilibrium conditions in real-experiments is often challenging. It often struggles for understanding phenomena like excited states or specific Arrhenius-driven reactions. As a result, the advantages of plasma modeling with simplifications can sometimes overshadow the disadvantages of experiments. This study simulated the ionization process of argon plasma using the 4th order Runge-Kutta numerical method. The simulation, initiated with initial densities before the simulation is run, each of them is electrons 2.6 × 1018 m-3, neutral argon (Ar) 2.6 × 1018 m-3, positive argon ions (Ar+) 2.6 × 1018 m-3, and positive diatomic argon ions (Ar2+) 2.6 × 1018 m-3, successfully obtained reaction rate equilibrium data at the 625th iteration. The final densities observed were 2.46 × 1018 m-3 for electrons, 2.27 × 1018 m-3 for neutral argon, 6.4 × 1015 m-3 for Ar+, and 4.34 × 1017 m-3 for Ar2+. These results show the equilibrium reaction rate in argon plasma which provides information that density of electron and Ar+ species show a decreasing trend while density of Ar and Ar2+ species shows an increasing trend which are the result of ionization and recombination processes in the entire plasma system.