We propose a numerical technique based on a combination of short-iterative Lanczos and exact diagonalization methods, suitable for simulating the time evolution of the reduced density matrix of a single qubit interacting with an environment. By choosing a mode discretization method and a flexible bath states truncation scheme, we are able to include in the physical description multiple-excitation processes, beyond weak coupling and Markov approximations. We apply our technique to the simulation of three different model Hamiltonians, which are relevant in the field of adiabatic quantum computation. We compare our results with those obtained on the basis of the widely used Lindblad master equation, as well as with well-known exact and approximated approaches. We show that our method is able to recover the thermodynamic behavior of the qubit-bath system, beyond the Born-Markov approximation. Finally, we show that even in the case of the adiabatic quantum annealing of a single qubit the bath can be beneficial in reaching the reduced system ground state.