ABSTRACT This study investigates the cosmic bulk flow through an analysis of luminosity distance variations in a perturbed Friedmann universe governed by f(R) gravity. The Hubble parameter, derived from perturbed Friedmann equations capturing intrinsic temporal fluctuations, is crucial. Redshift tomography is used to analyse observational data from the Pantheon catalogue from 0.015 to z < 2.3. Our objective is to constrain the cosmic bulk flow’s direction and magnitude within this redshift range. At low redshifts (z < 0.06), the predominant bulk flow aligns with the dominating supercluster in the corresponding range, maintaining a relatively constant magnitude, consistent with Lambda cold dark matter (ΛCDM) model predictions. A shift towards (l, b) = (290 ± 21, 15 ± 20) with $v_{\mathrm{bulk}} = 774 \pm 83 \, \mathrm{km\, s}^{-1}$ occurs between 0.06 < z < 0.1, coinciding with the cosmic microwave background dipole. For 0.1 < z < 0.2, the bulk flow shifts to (l, b) = (270 ± 21, 25 ± 20) with $v_{\mathrm{bulk}} = 903 \pm 102 \, \mathrm{km\, s}^{-1}$, supporting Kashlinsky et al.’s result. Beyond z > 0.2, the bulk flow’s magnitude exceeds $1000 \, \mathrm{km\, s}^{-1}$, deviating from ΛCDM model expectations. Our findings suggest that matter density oscillations under f(R) gravity exhibit ΛCDM model properties at low redshifts. However, the microwave background anisotropy shows noticeable fluctuations at higher redshifts, validating prior findings of strong cosmic bulk flows at these scales.