  We present a systematic numerical iteration approach to study the evolution
properties of the spin-boson systems, which works well in whole coupling
regime. This approach involves the evaluation of a set of coefficients for the
formal expansion of the time-dependent Schr\"{o}dinger equation $\vert
t\rangle=e^{-i\hat{H}t}\vert t=0\rangle$ by expanding the initial state $\vert
t=0\rangle$ in Fock space. The main advantage of this method is that this set
of coefficients is unique for the Hamiltonian studied, which allows one to
calculate the time evolution based on the different initial states. To
complement our numerical calculations, the method is applied to the
Buck-Sukumar model. Furthermore, we pointed out that, when the ground state
energy of the model is unbounded and no ground state exists in a certain
parameter space, the unstable time evolution of the physical quantities is the
natural results. Furthermore, we test the performance of the numerical method
to the Hamiltonian that use anti-Hermitian terms for modeling open quantum
systems.
