  We apply quadratic $f(R)=R+\epsilon R^2$ field equations, where $\epsilon$
has a dimension [L$^2$], to static spherical stellar model. We assume the
interior configuration is determined by Krori-Barua ansatz and additionally the
fluid is anisotropic. Using the astrophysical measurements of the pulsar PSR
J0740+6620 as inferred by NICER and XMM observations, we determine
$\epsilon\approx \pm 3$ km$^2$. We show that the model can provide a stable
configuration of the pulsar PSR J0740+6620 in both geometrical and physical
sectors. We show that the Krori-Barua ansatz within $f(R)$ quadratic gravity
provides semi-analytical relations between radial, $p_r$, and tangential,
$p_t$, pressures and density $\rho$ which can be expressed as $p_r\approx v_r^2
(\rho-\rho_1)$ and $p_r\approx v_t^2 (\rho-\rho_2)$, where $v_r$ ($v_t$) is the
sound speed in radial (tangential) direction, $\rho_1=\rho_s$ (surface density)
and $\rho_2$ are completely determined in terms of the model parameters. These
relations are in agreement with the best-fit equations of state as obtained in
the present study. We further put the upper limit on the compactness, which
satisfies the $f(R)$ modified Buchdahl limit. Interestingly, the quadratic
$f(R)$ gravity with negative $\epsilon$ naturally restricts the maximum
compactness to values lower than Buchdahl limit, unlike the GR or $f(R)$
gravity with positive $\epsilon$ where the compactness can arbitrarily approach
the black hole limit $C\to 1$. The model predicts a core density a few times
the saturation nuclear density $\rho_{\text{nuc}} = 2.7\times 10^{14}$
g/cm$^3$, and a surface density $\rho_s > \rho_{\text{nuc}}$. We provide the
mass-radius diagram corresponding to the obtained boundary density which has
been shown to be in agreement with other observations.
