14 (102) or (103) or (105) can be handled in practice as constrained ordinary differential equations, solving nu- merically (106) can be achieved with simple numerical techniques. Since δ ( H ) diverges as H q → 0, as the maximal expan- sion ct = ct max is reached in elliptic models or regions (see A 4), it is better to use in these cases the alterna- tive variable Σ = −H $_{q}$δ ( H ) (see (32)). This leads to the alternative system We may supplement either (102) or (103) with the fol- lowing extra pair of differential equations where we used (31) and R$_{i}$ can be prescribed as a radial coordinate gauge. An obvious choice of initial value functions for the sys- tems (102)–(104) or (103)–(104) is B. Initial conditions. m$_{qi}$, H $_{qi}$, δ ( m ) i = m ′ $_{qi}$/m$_{qi}$ 3 R ′ $_{i}$/R$_{i}$ , δ ( H ) i = H ′ $_{qi}$/ H qi 3 R ′ $_{i}$/R$_{i}$ or Σ$_{i}$ = − H ′ qi 3 R ′ $_{i}$/R$_{i}$ (107) ˙ m$_{q}$ = − 3 m$_{q}$ H $_{q}$, (103a) ˙ H q = −H 2 q − m$_{q}$, (103b) ˙ δ ( m ) = − 3(1 + δ ( m $^{)}$) Σ , (103c) ˙ Σ = 2 Σ H q + Σ 2 + m$_{q}$ δ ( m $^{)}$, (103d) which is entirely equivalent to (102), and has the advan- tage (besides its use for recollapsing configurations) that Σ is the scalar associated with the shear tensor by (15). which follow directly from (32), (38) and (40). By includ- ing these extra equations we can use the systems above to determine numerically the metric functions in (39). In fact, if we substitute the scaling laws (37)–(38), (41), (42) and (43) into (104) we can work only with these two evolution equations given as Since k$_{qi}$ determines the solutions of (38) (or (3) for R = R$_{i}$L ), this initial value function is closely related to the kinematic evolution of the models (or specific regions of them). Hence, it is more practical and intuitive to choose instead of (107) the initial value functions ˙ L = L H $_{q}$, (104a) ˙ Γ = 3 Γ $_{q}$δ ( H ) = − 3 Γ Σ , (104b) m$_{qi}$, k$_{qi}$, δ ( m ) i = m ′ $_{qi}$/m$_{qi}$ 3 R ′ $_{i}$/R$_{i}$ , δ ( k ) i = k ′ $_{qi}$/k$_{qi}$ 3 R ′ $_{i}$/R$_{i}$ . (108) These initial conditions are sufficient to solve (105), but to solve (102)–(104) or (104)–(103) we also need H qi and δ ( H ) i or Σ$_{i}$, which readily follow from (38) and (43) whose solution ( L and Γ) allows us to compute all scalars by means of these scaling laws. However, from a numeric point of view it might be more convenient and practical to work with (102) or (103) supplemented by (104) than with (105). Solving any one of the systems of evolution equations described above is needed to find radial null geodesics and their corresponding red shift factor, z , which are in turn needed for computing observable quantities (lumi- nosity distance, etc). Given a numerical solution of either system, radial geodesics ct = ct$_{N}$ ( r ) in the past null cone and z follow from solving [3] The main advantage of using (108) as initial conditions is the fact that they are exactly the same initial value functions that we employed in the parametrization of the analytic solutions and in the forms of Γ in (59)–(60), and the Hellaby–Lake conditions. Hence, given a choice of these functions, a time evolution free from shell crossings can be immediately tested by means of (64), (68) and (72). The evolution equations (103) and have already been used for a dynamical systems approach to LTB mod- els [29], while a suitable generalization of (102) was employed for a numeric study of cosmological models endowed with the LTB metric but with an energy– momentum tensor of an anisotropic fluid [31]. A comprehensive numeric study of LTB models by means of the systems (102), (103) or (105) is beyond the scope of this article, and thus is presently under ex- amination in a separate work. Our purpose has been to show how these systems, constructed with quasi–local variables, are potentially promising for this purpose. H qi = [2 m$_{qi}$ − k$_{qi}$ ] 1 / $^{2}$, δ ( H ) i = 2 m$_{qi}$δ ( m ) i − k$_{qi}$δ ( k ) i 2[2 m$_{qi}$ k$_{qi}$ ] or Σ$_{i}$ = − 2 m$_{qi}$δ ( m ) i − k$_{qi}$δ ( k ) i 2[2 m$_{qi}$ k$_{qi}$ ] 1 / 2 , (109) c d t$_{N}$ ( r ) d r = − R ′ N (1 + E ) 1 / 2 = − L$_{N}$ Γ$_{N}$ R ′ i [1 − k$_{qi}$R 2 $_{i}$] 1 / $_{2}$, (106a) d d r ln(1 + z ) = ˙ R ′ N (1 + E ) 1 / 2 = L$_{N}$ Γ$_{N}$ H qN (1 + 3 δ ( H ) N ) R ′ i [1 − k$_{qi}$R 2 $_{i}$] 1 / 2 , (106b) where the subindex N means evaluation at the null curve ct = ct$_{N}$ ( r ), i.e. A$_{N}$ = A ( ct$_{N}$ ( r ) , r ). Since the systems