  The first part of this paper introduces sufficient conditions to determine
conservation laws of diffusion equations of arbitrary fractional order in time.
Numerical methods that satisfy a discrete analogue of these conditions have
conservation laws that approximate the continuous ones. In the second part of
the paper, we propose a method that combines a finite difference method in
space with a spectral integrator in time. The time integrator has already been
applied in literature to solve time fractional equations with Caputo fractional
derivative of order $\alpha\in(0,1)$. It is here generalised to approximate
Caputo and Riemann-Liouville fractional derivatives of arbitrary order. We
apply the method to subdiffusion and superdiffusion equations with
Riemann-Liouville fractional derivative and derive its conservation laws.
Finally, we present a range of numerical experiments to show the convergence of
the method and its conservation properties.
