I. INTRODUCTION It is one of the most important issues for cosmologists and particle physicists to under- stand the physical origin of the dark energy (DE) which is responsible for an accelerated expansion of the current Universe. Although the standard spatially flat Λ-Cold-Dark-Matter (ΛCDM) model is consistent with all kinds of current observational data [1], some tentative deviations from it have been reported recently [2, 3]. Furthermore, in the ΛCDM model, the cosmological term is regarded as a new fundamental constant whose observed value is much smaller than any other energy scale known in physics. Hence it is natural to seek for non-stationary models of the current DE. Among them, f ( R ) gravity which modifies and generalizes the Einstein gravity by incorporating a new phenomenological function of the Ricci scalar R , f ( R ), can provide a self-consistent and non-trivial alternative to the ΛCDM model[4, 5]. In the previous paper [6], we calculated evolution of matter density fluctuations in viable f ( R ) models [4, 5] for redshifts z ≫ 1 during the matter-dominated stage and found an analytic expression for them. In this paper we extend the previous analysis and perform numerical calculations of the evolution of both background space-time and density fluctu- ations for the particular f ( R ) model of Ref. [5] without such a restriction. As a result, we have found crossing of the phantom boundary w$_{DE}$ = − 1 at an intermediate redshift z ≲ 1 for the background space-time metric and an anomalous behavior of the growth index of fluctuations. We adopt the following action of f ( R ) models with model parameters n, λ and R$_{s}$ [5]: where S$_{m}$ is the action of the matter content which is assumed to be minimally coupled to gravity. To make the late-time asymptotic de Sitter regime where R = constant stable, λ has to satisfy f $^{′}$( R ) > Rf $^{′′}$( R ). As a result, λ has a lower limit λ$_{min}$ for each n . Numerically we find ( n, λ$_{min}$ ) =(2, 0.9440), (3, 0.7259), and (4, 0.6081). From the action (1), we obtain 2 S = 1 16 πG ∫ d 4 $_{x}$√ − gf ( R ) + S$_{m}$, f ( R ) = R + λR$_{s}$ [ ( 1 + R 2 R 2 s ) − n − 1 ] , (1) II. BACKGROUND