Dynamics of the sub-Ohmic spin-boson model is examined using three numerical approaches, namely the Dirac-Frenkel time-dependent variation with the Davydov D1 ansatz, the adaptive time-dependent density matrix renormalization group method within the representation of orthogonal polynomials, and a perturbative approach based on a unitary transformation. In order to probe the validity regimes of the three approaches, we study the dynamics of a qubit coupled to a bosonic bath with and without a local field. Comparison of the up-state population evolution shows that the three approaches are in agreement in the weak-coupling regime but exhibit marked differences when the coupling strength is large. The Davydov D1 ansatz and the time-dependent density matrix renormalization group can both be reliably employed in the weak-coupling regime, while the former is also valid in the strong-coupling regime as judged by how faithfully the trial state follows the Schrodinger ¨ equation. We further explore the bipartite entanglement dynamics between two qubits coupled with individual bosonic baths which reveals entanglement sudden death and revival.