We　investigate　the　suitability　of　a　fractional　Laplacian　in　the　construction　of　a　three-dimensional　(3D)　multi-term　time–space　fractional　Bloch–Torrey　equation.　Instead　of　using　the　Riesz　fractional　derivative　as　in　most　previous　research,　we　consider　the　fractional　Laplacian　in　space　and　multi-term　Caputo　fractional　derivative　in　time.　A　vertex-centred　finite　volume　method　is　firstly　derived　to　discretise　this　3D　fractional　dynamic　system　in　combination　with　the　matrix　transfer　technique　and　the　weighted　and　shifted　Grünwald–Letnikov　formula.　Since　computing　the　fractional　power　of　the　matrix　representation　m(−Δ)　of　the　Laplacian　operator　results　in　a　highly　dense　matrix,　we　use　matrix　function　approximations　developed　in　a　Krylov　subspace　framework　and　exploit　the　sparsity　of　m(−Δ)　to　reduce　the　computational　overhead.　Then　the　Lanczos　method　is　employed　to　approximate　vector　products　of　matrix　functions.　A　preconditioner　is　applied　to　accelerate　the　convergence　process,　which　shows　excellent　performance　and　reduces　the　memory　considerably.　In　particular,　we　introduce　block　matrix　functions　to　deal　with　the　coupled　system.　Furthermore,　the　analytical　solution　of　the　two-term　time–space　fractional　diffusion　equation　is　derived　and　expressed　by　the　multinomial　Mittag-Leffler　function.　The　feasibility　and　effectiveness　of　the　proposed　numerical　techniques　are　verified　by　some　numerical　examples.　Compared　with　two-dimensional　fractional　models　or　single-term　time　fractional　dynamic　systems,　this　multi-term　time–space　fractional　model　provides　more　flexibility　for　the　characterisation　of　anomalous　diffusion　in　biological　tissues.　©　2022　Elsevier　B.V.