$\def\ii{{\rm i}}
\def\dd{{\rm d}}
\def\ee{{\rm e}}
$
It turns out that the case of pure dephasing is exactly solvable, and one can obtain nice solutions under certain conditions. In particular, I will consider the case of Gaussian, stationary noise.
Exact solution
Let us define the noisy qubit Hamiltonian $(\hbar = 1)$
$$ \hat{H}(t) = \frac{1}{2}\left[\omega_0 + \delta \hat{\omega}(t) \right]\sigma^z, $$
where $\omega_0$ represents the fixed qubit resonance frequency and $\delta\hat{\omega}(t)$ is a zero-mean random variable describing the effect of some external fluctuating fields. (Hats indicate random variables). Because the Hamiltonian commutes with itself at all times, the time evolution operator can be calculated exactly:
$$ \hat{U}(t) = T\exp \left( -\ii \int_0^t\dd s\; \hat{H}(s) \right) = \exp \left( -\frac{\ii}{2}\left[\omega_0 t +\hat{\Phi}(t) \right]\sigma^z\right),$$
where we defined a new random variable $$\hat{\Phi}(t) = \int_0^t\dd s\;\delta\hat{\omega}(s). $$
The solution for the density matrix given some fixed initial condition $\rho(0)$ is of course $\hat{\rho}(t) = \hat{U}(t)\rho(0) \hat{U}^\dagger(t)$. The populations of $\hat{\rho}(t)$ in the eigenbasis of $\sigma^z$ are constant, while the off-diagonal elements $\hat{\rho}_{01}(t) = \hat{\rho}_{10}^*(t)$ are given by $$\hat{\rho}_{01}(t) = \ee^{\ii\omega_0 t}\ee^{\ii\hat{\Phi}(t)} \rho_{01}(0). $$
Given a realisation of the stochastic process $\delta\hat{\omega}(t)$, one can compute the corresponding coherence $\hat{\rho}_{01}(t)$ using the above. Note that $\hat{\rho}_{01}$ is also a random variable describing a typical evolution of a single qubit. More useful, however, is a description of the ensemble obtained by averaging over many identically prepared copies of the system, since this is what one actually does in the experiment. Thus, what we really want is
$$\left \langle \ee^{\ii\hat{\Phi}(t)} \right\rangle = \ee^{-\Gamma(t)},$$
where the angle brackets denote an average over many realisations of the noise, and the formal object $\Gamma(t)$ is called the decoherence function. It follows from the above that the ensemble-averaged density matrix $\rho(t) = \langle \hat{\rho}(t)\rangle$ obeys the exact, non-Markovian master equation
$$ \dot{\rho} = -\frac{\ii}{2}\left(\omega_0 - \Sigma(t)\right)[\sigma^z,\rho] + \frac{1}{2}\gamma(t)\left( \sigma^z \rho \sigma^z - \rho \right), $$
where the dot denotes a time derivative and we defined the energy shift $\Sigma(t) = {{\rm Im}}\left[\dot{\Gamma}\right]$ and the decoherence rate $\gamma(t) = {{\rm Re}}\left[\dot{\Gamma}\right]$.
Gaussian, stationary noise
In order to make progress, we need to make some assumptions about the noise. First, we recognise that the decoherence function is (up to a sign) equal to the cumulant-generating function of $\hat{\Phi}(t)$. Now, if we assume that $\hat{\Phi}(t)$ is Gaussian, then only its second cumulant is non-zero (the first cumulant vanishes because $\langle \delta\hat{\omega}(t)\rangle = 0$). It follows that the decoherence function is given by the autocorrelation function of the fluctuating frequency:
$$ \Gamma(t) = \frac{1}{2}\int_0^t\dd s\int_0^t\dd s'\; \langle \delta\hat{\omega}(s) \delta\hat{\omega}(s')\rangle.$$
Let us further assume that the noise statistics are stationary, i.e. invariant under time translations, which implies that $\langle \delta\hat{\omega}(s) \delta\hat{\omega}(s')\rangle = f(s-s')$, and therefore
$$\Gamma(t) = \frac{1}{2} \int_0^t \dd\tau\int_{-\tau}^\tau\dd\tau' f(\tau'),$$
where we changed variables to $\tau = s+s'$ and $\tau' = s-s'$. Clearly, $\Sigma(t) = 0$, while a few simple manipulations should convince you that
\begin{align} \gamma(t) & = \frac{1}{2} \int_{-t}^t \dd\tau\;f(\tau) \qquad \qquad (*)\\
& = \int_{-\infty}^{\infty} \dd \omega\; \frac{\sin(\omega t)}{\omega} S(\omega),
\end{align}
where we defined the noise power spectral density
$$ S(\omega) = \frac{1}{2\pi}\int_{-\infty}^{\infty} \dd \tau\; \ee^{-\ii \omega \tau} f(\tau).$$
Now, for almost any realistic noise source, the autocorrelation function $f(\tau)$ must eventually decay to zero when $|\tau|$ is greater than some time scale denoted $T$. For times longer than $T$, one may extend the integration limits to $\pm \infty$ in the equation marked $(*)$ above. This leads to an asymptotic decay rate $\gamma_\infty = \lim_{t\to \infty} \gamma(t) = \pi S(0)$.
The Markov approximation consists of the assumption that this time scale $T$ is much shorter than the typical time scale of the qubit's noisy evolution, which here is fixed by $\gamma(t)$ (or rather its inverse). This is equivalent to saying that the noise power spectral density varies slowly over a frequency scale $\sim\gamma(t)$. In this case, it is possible to immediately replace $\gamma(t)\to \gamma_\infty$ with negligible error, leading to simple exponential decoherence $\rho_{01}(t) \sim\ee^{-\gamma_\infty t}$. This is most easily illustrated with a simple example, as follows.
Example 1
Let us take a simple exponential form for the noise autocorrelation function
$$\langle \delta\hat{\omega}(s) \delta\hat{\omega}(s')\rangle = g^2\ee^{-\Omega|s-s'|},$$
This clearly falls to zero over a time scale $T =1/\Omega$. The decoherence rate is
$$ \gamma(t) = \frac{g^2}{\Omega} ( 1- \ee^{-\Omega t}),$$
and the noise power spectral density is
$$ S(\omega) = \frac{1}{\pi} \frac{g^2\Omega}{\omega^2 + \Omega^2},$$
implying an asymptotic decoherence rate $\gamma_\infty = \pi S(0) = g^2/\Omega$. The dynamics will be approximately Markovian if $\Omega \gg \gamma_\infty$, or in other words if $\Omega\gg g$. This means that the coupling of the qubit to its noisy environment ($g$) must be small in comparison to the bandwidth of the noise ($\Omega$).
The plots below show the exact evolution of (the absolute value of) the qubit coherence for two examples, with the Markov approximation included for comparison (dotted line). The first plot shows $\Omega = 2\gamma_\infty$, where the transient dynamics differs significantly from the simple Markov result. The second plot shows $\Omega = 100\gamma_\infty$, in which case the noise memory time is negligibly short and the Markov approximation is accurate for all practical purposes. Note, however, that there always exists a time regime $t\ll \Omega^{-1}$ where the dynamics is non-Markovian. Nevertheless, whenever the Markov assumption $\gamma_\infty\ll \Omega$ holds, the system barely evolves over such a short time and thus the non-Markovian aspects of the evolution are not seen.
Example 2
Just for fun, here's a more complex example. Consider the following spectral density
$$ S(\omega) = \frac{\alpha}{\Omega^2} \omega^3 \left[ 1 - \mathrm{sinc}\left( \frac{\omega}{\varpi}\right) \right] \ee^{-\omega^2/(2\Omega^2)}.$$
Such a spectral density can be realised with a qubit comprising two localised states of an atom in a double-well potential. If this atom is plunged into a cold bosonic gas forming a Bose-Einstein condensate, it decoheres by scattering atoms in the condensate and thus generating phonons. In particular, if the impurity atom (i.e. the one forming the qubit) is prepared in a delocalised state, i.e. a coherent superposition of being localised in the left well ($\lvert 0\rangle$) and the right well ($\lvert 1\rangle$), over time this coherence vanishes. The coupling $\alpha$ is proportional to the scattering cross-section for impurity-gas collisions, while $\Omega = c/\ell$ and $\varpi = c/L$, where $\ell$ is the width of the wave functions describing the localised states, $L$ is the separation between the double-well minima, and $c$ is the speed of sound in the gas. More details can be found in, for example, one of our papers here.
See below for a plot of the decoherence dynamics, for $\alpha = 1$ and $\Omega = 10\varpi$. Interestingly, the coherence does not asymptotically vanish, but reaches a constant, finite value. Even more interestingly, there are two periods where lost coherence is regained, i.e. the decoherence rate is negative, a signature of truly non-Markovian dynamics.
The second recoherence phase (occurring at $\Omega t = 10$) has a particularly pleasing physical interpretation: it always occurs at a time $t = \varpi^{-1} = L/c$, which is the time taken for the disturbance in the superfluid generated at one potential minimum to propagate to the other one. Thus some of the lost coherence is regained due to the fact that both localised states interact with a common superfluid reservoir. However, this only comes into effect after some time delay due to the finite speed of signal propagation in the superfluid.
