  In this article, two kinds of numerical algorithms are derived for the
ultra-slow (or superslow) diffusion equation in one and two space dimensions,
where the ultra-slow diffusion is characterized by the Caputo-Hadamard
fractional derivative of order $\alpha \in (0,1)$. To describe the spatial
interaction, the Riesz fractional derivative and the fractional Laplacian are
used in one and two space dimensions, respectively. The Caputo-Hadamard
derivative is discretized by two typical approximate formulae, i.e.,
L2-1$_{\sigma}$ and L1-2 methods. The spatial fractional derivatives are
discretized by the 2-nd order finite difference methods. When L2-1$_{\sigma}$
discretization is used, the derived numerical scheme is unconditionally stable
with error estimate $\mathcal{O}(\tau^{2}+h^{2})$ for all $\alpha \in (0, 1)$,
in which $\tau$ and $h$ are temporal and spatial stepsizes, respectively. When
L1-2 discretization is used, the derived numerical scheme is stable with error
estimate $\mathcal{O}(\tau^{3-\alpha}+h^{2})$ for $\alpha \in (0, 0.3738)$. The
illustrative examples displayed are in line with the theoretical analysis.
