  For cosmologically interesting $f(R)$ gravity models, we derive the complete
set of the linearized field equations in the Newtonian gauge, under
environments of the solar system, galaxies and clusters respectively. Based on
these equations, we confirmed previous $\gamma=1/2$ solution in the solar
system. However, $f(R)$ gravity models can be strongly environment-dependent
and the high density (comparing to the cosmological mean) solar system
environment can excite a viable $\gamma=1$ solution for some $f(R)$ gravity
models. Although for $f(R)\propto -1/R$, it is not the case; for $f(R)\propto
-\exp(-R/\lambda_2H_0^2)$, such $\gamma=1$ solution does exist. This solution
is virtually indistinguishable from that in general relativity (GR) and the
value of the associated curvature approaches the GR limit, which is much higher
than value in the $\gamma=1/2$ solution. We show that for some forms of $f(R)$
gravity, this solution is physically stable in the solar system and can
smoothly connect to the surface of the Sun.
  The derived field equations can be applied directly to gravitational lensing
of galaxies and clusters. We find that, despite significant difference in the
environments of galaxies and clusters comparing to that of the solar system,
gravitational lensing of galaxies and clusters can be virtually identical to
that in GR, for some forms of $f(R)$ gravity. Fortunately, galaxy rotation
curve and intra-cluster gas pressure profile may contain valuable information
to distinguish these $f(R)$ gravity models from GR.
