The point kinetics equations with two-energy groups and multi-group of delayed neutrons are a system of stiff linear/nonlinear ordinary differential equations strongly coupled, which is used to simulate the dynamics system of a homogeneous nuclear reactor by determining the population of neutrons and precursor groups. The matrix form of the mathematical model is solved by the theta method where the matrix inverse is calculated analytically. Numerical results are obtained through iterative schemes to approximate analytical solutions. The stiffness condition of the system of equations can cause instability in the solutions of explicit iterative schemes commonly employed so far. In this work, an implicit iterative scheme is used to approximate numerical solutions to the homogeneous nuclear reactor with one and six groups of delayed neutrons. The results are discussed and compared on light of literature reports. Unprecedented comparisons with experimental data from an AGN-201 nuclear reactor are also discussed. For the case studies, the theta method exhibits accurate and computationally efficient results.