We evaluate the non-Markovian finite-temperature two-time correlation functions (CFs) of system operators of a pure-dephasing spin-boson model in two different ways, one by the direct exact operator technique and the other by the recently derived evolution equations, valid to second order in the system-environment interaction Hamiltonian. This pure-dephasing spin-boson model that is exactly solvable has been extensively studied as a simple decoherence model. However, its exact non-Markovian finite-temperature two-time system operator CFs, to our knowledge, have not been presented in the literature. This may be mainly due to the fact, illustrated in this article, that in contrast to the Markovian case, the time evolution of the reduced density matrix of the system (or the reduced quantum master equation) alone is not sufficient to calculate the two-time system operator CFs of non-Markovian open systems. The two-time CFs obtained using the recently derived evolution equations in the weak system-environment coupling case for this non-Markovian pure-dephasing model happen to be the same as those obtained from the exact evaluation. However, these results significantly differ from the non-Markovian two-time CFs obtained by wrongly directly applying the quantum regression theorem (QRT), a useful procedure to calculate the two-time CFs for weak-coupling Markovian open systems. This demonstrates clearly that the recently derived evolution equations generalize correctly the QRT to non-Markovian finite-temperature cases. It is believed that these evolution equations will have applications in many different branches of physics.