results are presented in section IV. We give our conclusions and discussions in section V. To fit models with observational data, we need to know the form of Hubble function H ( z ) (or E ( z ) = H ( z ) /H$_{0}$ ). In a flat FRW universe where Ω (0) r , Ω (0) b , Ω (0) dm and Ω (0) de are present values of the dimensionless energy density for radiations, baryons, dark matter and dark energy, respectively, and Ω (0) r + Ω (0) b + Ω (0) dm + Ω (0) de = 1. The energy densities of baryons and dark matter are always written together as Ω (0) b (1 + z ) 3 + Ω (0) $_{dm}$(1 + z ) 3 = Ω (0) m (1 + z ) $^{3}$. The radiation density is the sum of photons and relativistic neutrinos [20]: where N$_{n}$ is the number of neutrino species and Ω (0) γ = 2 . 469 × 10 − $^{5}$h − 2 for T$_{cmb}$ = 2 . 725 K ( h = H$_{0}$/ 100 Mpc · km · s − $^{1}$). The evolving function F ( z ) for DE depends on w$_{de}$ ( z ): For example, for the CPL parametrization and F ( z ) = 1 for w$_{de}$ = − 1, respectively. Here we divide z ∈ (0 , ∞ ) into m + 1 bins and assume w$_{de}$ ( z ) in the first m bins as and require w$_{de}$ ( z ) to be continuous at divided points: Note that here prime does not represent a derivative, instead w ′ n is just the slope of the linear expansion in the n th bin. Thus the independent parameters are w$_{0}$, w ′ $_{1}$, w ′ $_{2}$, ..., w ′ $_{n}$, ..., w ′ m (5) w$_{n}$ = w$_{n}$$_{−}$$_{1}$ + w ′ n × ( z$_{n}$ − z$_{n}$$_{−}$$_{1}$ ) , (1 ≤ n ≤ m − 1) (4) w$_{de}$ ( z$_{n}$$_{−}$$_{1}$ < z ≤ z$_{n}$ ) = w$_{n}$$_{−}$$_{1}$ + w ′ n × ( z − z$_{n}$$_{−}$$_{1}$ ) , (1 ≤ n ≤ m ) (3) F ( z ) = (1 + z ) 3(1+ w$_{0}$ + w a $^{)}$e − 3 w a z 1+ z , F ( z ) = e 3 ∫ z 0 1+ w de 1+ x $^{dx}$. (2) Ω (0) r = Ω (0) γ (1 + 0 . 2271 N$_{n}$ ) , E $^{2}$( z ) = Ω (0) r (1 + z ) 4 + Ω (0) b (1 + z ) 3 + Ω (0) $_{dm}$(1 + z ) 3 + Ω (0) $_{de}$F ( z ) , (1) II. METHODOLOGY 4