Magnetic resonance imaging can reveal exquisite details about the complex structure and function of human tissue. However, magnetic resonance imaging signal behaviour at high or ultra-high field has shown increased deviation from the classically expected mono-exponential relaxation. The underlying mechanism of anomalous relaxation can contribute to a better understanding of the interaction of spins with their surroundings. The purpose of this work is to explore the utility of the multi-term time-fractional Bloch equations in relation to the anomalous relaxation processes. We proposed an effective predictor–corrector method to solve the multi-term equations. Voxel-level temporal fitting of the magnetic resonance imaging signal was performed based on the model developed from the multi-term time-fractional Bloch equations. A feasible parameter estimation method based on hybrid Nelder–Mead simplex search and particle swarm optimisation was introduced to perform the curve fitting. The extra time-fractional terms provide flexibility in the relaxation process.