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 α ∈ (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.