We consider the Bloch–Torrey operator, − Δ + i g x , that governs the time evolution of the transverse magnetization in diffusion magnetic resonance imaging (dMRI). Using the matrix formalism, we compute numerically the eigenvalues and eigenfunctions of this non-Hermitian operator for two bounded three-dimensional domains: a sphere and a capped cylinder. We study the dependence of its eigenvalues and eigenfunctions on the parameter g and on the shape of the domain (its eventual symmetries and anisotropy). In particular, we show how an eigenfunction drastically changes its shape when the associated eigenvalue crosses a branch (or exceptional) point in the spectrum. Potential implications of this behavior for dMRI are discussed.