Fractional differential equations play an important role in science and technology. Many problems can be cast using both fractional time and spatial derivatives. In order to accurately simulate natural phenomena using this technology one needs fine spatial and temporal discretizations. This leads to large-scale linear systems or matrix equations, especially whenever more than one space dimension is considered. The discretization of fractional differential equations typically involves dense matrices with a Toeplitz structure. We combine the fast evaluation of Toeplitz matrices and their circulant preconditioners with state-of-the-art linear matrix equation solvers to efficiently solve these problems, both in terms of CPU time and memory requirements. Numerical experiments on typical differential problems with fractional derivatives in both space and time showing the effectiveness of the approaches are reported.