We present an application of the extended stochastic Liouville-von Neumann equation (ESLN) method introduced earlier [G. M. G. McCaul, C. D. Lorenz, and L. Kantorovich, Phys. Rev. B 95, 125124 (2017)2469-995010.1103/PhysRevB.95.125124; G. M. G. McCaul, C. D. Lorenz, and L. Kantorovich, Phys. Rev. B 97, 224310 (2018)]2469-995010.1103/PhysRevB.97.224310, which describes the dynamics of an exactly thermalized open quantum system reduced density matrix coupled to a non-Markovian harmonic environment. Critically, the combined system of the open system fully coupled to its environment is thermalized at finite temperature using an imaginary-time evolution procedure before the application of real-time evolution. This initializes the combined system in the correct canonical equilibrium state rather than being initially decoupled. Here we apply our theory to the spin-boson Hamiltonian and develop a number of competing ESLN variants designed to reduce the numerical divergence of the trace of the open-system density matrix. We find that a careful choice of the driving noises is essential for improving numerical stability. We also investigate the effect of applying higher-order numerical schemes for solving stochastic differential equations, such as the Stratonovich-Heun scheme, and conclude that stochastic sampling dominates convergence with the improvement associated with the numerical scheme being less important for short times but required for late times. To verify the method and its numerical implementation, we first consider evolution under a fixed Hamiltonian and show that the system either remains in, or approaches, the correct canonical equilibrium state at long times. Additionally, evolution of the open system under nonequilibrium Landau-Zener (LZ) driving is considered and the asymptotic convergence to the LZ limit was observed for vanishing system-environment coupling and temperature. When coupling and temperature are nonzero, initially thermalizing the combined system at a finite time in the past was found to be a better approximation of the true LZ initial state than starting in a pure state.