  Modification of general relativity (GR) inspired by theories like $f(R)$
gravity is among the most popular ones to explain the late-time acceleration of
the Universe as an alternative to the $\Lambda$CDM model. In this work, we use
the state-of-the-art BAO+BBN data and the most recent Type Ia supernovae (SNe
Ia) sample namely PantheonPlus, including the Cepheid host distances and
covariance from SH0ES samples, to robustly constrain the $f(R)$ gravity
framework via two of the most popular $f(R)$ models in literature, namely, the
Hu-Sawicki and Starobinsky models. Additionally, we consider how the time
variation of the Newton's gravitational constant affects the supernovae
distance modulus relation. We find a minor evidence for $f(R)$ gravity under
the Hu-Sawicki dynamics from BAO+BBN and BAO+BBN+uncalibrated supernovae joint
analysis, but the inclusion of Cepheid host distances, makes the model
compatible with GR. Further, we notice tendency of this model to relax the
$H_0$ tension. In general, in all the analyses carried out in this study with
the late time probes, we find both the $f(R)$ models to be consistent with GR
at 95\% CL.

Cosmology and Nongalactic Astrophysics;General Relativity and Quantum Cosmology