  In the present paper, we treat the problem of the existence of quark stars
(QSs) for selected homogeneous and unpaired charge-neutral $3$-flavor
interacting quark matter with $\mathcal{O}(m_s^4)$ corrections equations of
state (EoS). Using the EoS combined with the Tolmann-Oppenheimer-Volkoff (TOV)
structure equations, the properties of stars are explored by obtaining their
mass-radius relations. All calculations are carried out within the framework of
the $R$-squared gravity defined by $f(R) = R+a R^2$. Our main goal is to
discuss the effect for a wide range of the $R$-squared gravity parameter, $a$,
on the mass-radius and the mass-central mass density $(M-\epsilon_c)$ relation
of QSs. Furthermore, we investigate the dynamical stability condition for those
stars, and we show that their dynamical stability depends on the anisotropic
parameters $ B_{\perp}$ and $a_4^{\perp}$ coming from the respective EoS. In
such a scenario our results provide circumstantial evidence in favor of
super-massive pulsars.
