  We perform a maximum likelihood kinematic analysis of the two dynamically
relaxed galaxy clusters MACS J1206.2-0847 at $z=0.44$ and RXC J2248.7-4431 at
$z=0.35$ to determine the total mass profile in modified gravity models, using
a modified version of the MAMPOSSt code of Mamon, Biviano and Bou\'e. Our work
is based on the kinematic and lensing mass profiles derived using the data from
the Cluster Lensing And Supernova survey with Hubble (hereafter CLASH) and the
spectroscopic follow-up with the Very Large Telescope (hereafter CLASH-VLT). We
assume a spherical Navarro-Frenk-White (NFW hereafter) profile in order to
obtain a constraint on the fifth force interaction range $\lambda$ for models
in which the dependence of this parameter on the enviroment is negligible at
the scale considered (i.e. $\lambda=const$) and fixing the fifth force strength
to the value predicted in $f(R)$ gravity. We then use information from lensing
analysis to put a prior on the other NFW free parameters. In the case of MACSJ
1206 the joint kinematic+lensing analysis leads to an upper limit on the
effective interaction range $\lambda\le1.61\,\mbox{Mpc}$ at
$\Delta\chi^{2}=2.71$ on the marginalized distribution. For RXJ 2248 instead a
possible tension with the $\Lambda$CDM model appears when adding lensing
information, with a lower limit $\lambda\ge0.14\,\mbox{Mpc}$ at
$\Delta\chi^{2}=2.71$. This is consequence of the slight difference between the
lensing and kinematic data, appearing in GR for this cluster, that could be
explained in terms of modifications of gravity. We discuss the impact of
systematics and the limits of our analysis as well as future improvements of
the results obtained.[...]
