FIG. 1: Numerical (solid) and approximate (dot-dashed) solutions for (a) all D 3 ( A ), D 2 ( A ), and D 1 ( A ), and (b) the enlarged plot of D 2 ( A ) and D 1 ( A ) in 0 ≤ A/δ 0 ≤ 10. 10 20 30 40 50 60 A GLYPH<Slash1> Δ$_{0}$ 10 20 30 40 50 60 D$_{i}$ GLYPH<Slash1> Δ$_{0}$ D$_{3}$ D$_{2}$ D$_{1}$ GLYPH<LParen1> a GLYPH<RParen1> 2 4 6 8 10 A GLYPH<Slash1> Δ$_{0}$ 2 4 6 8 10 D$_{i}$ GLYPH<Slash1> Δ$_{0}$ D$_{1}$ D$_{2}$ GLYPH<LParen1> b GLYPH<RParen1> 10 GLYPH<Minus> 3 10 GLYPH<Minus> 2 10 GLYPH<Minus> 1 1 10 10 2 10 $_{3}$A GLYPH<Slash1> Δ$_{0}$ 0.05 0.1 0.15 0.2 0.25 0.3 x$_{1}$ 10 GLYPH<Minus> 3 10 GLYPH<Minus> 2 10 GLYPH<Minus> 1 1 10 10 2 10 $_{3}$A GLYPH<Slash1> Δ$_{0}$ 0.1 0.2 0.3 0.4 x$_{2}$ 10 GLYPH<Minus> 3 10 GLYPH<Minus> 2 10 GLYPH<Minus> 1 1 10 10 2 10 $_{3}$A GLYPH<Slash1> Δ$_{0}$ 0.1 0.2 0.3 0.4 0.5 x$_{3}$ FIG. 2: The numerical (solid) and approximate (dot-dashed) solutions for x 1 ( A ), x 2 ( A ), and x 3 ( A ). Note that y i ( A ) ≃ − x i ( A ). W$_{0}$ ∼ =    2 / 3 1 / 3 0 1 / 6 1 / 3 1 / 2 1 / 6 1 / 3 1 / 2    , W$_{l}$ ∼ =    1 / 2 1 / 2 0 1 / 4 1 / 4 1 / 2 1 / 4 1 / 4 1 / 2    , W$_{i}$ ∼ =    0 1 0 1 / 2 0 1 / 2 1 / 2 0 1 / 2    , W$_{h}$ ∼ =    0 1 / 2 1 / 2 1 / 2 1 / 4 1 / 4 1 / 2 1 / 4 1 / 4    , W$_{d}$ ∼ =    0 0 1 1 / 2 1 / 2 0 1 / 2 1 / 2 0    . (28) A = ( A$_{0}$, A$_{l}$, A$_{i}$, A$_{h}$, A$_{d}$ ): Together, these matrices exhibit the remarkable simplic- ity of the PMNS matrix as A varies from 0 to ∞ . Note that all of the matrices have at least one zero, W$_{1}$$_{I}$ = 0, implying x$_{I}$ = y$_{I}$ = 0. Also, they have equal elements in their second and third rows, W$_{2}$$_{i}$ = W$_{3}$$_{i}$ , so that w$_{1}$$_{i}$ = 0 or x$_{i}$ + y$_{i}$ = 0. As a consequence, using the unitarity conditions, the W matrix is completely fixed by its el- ements in the first row, W$_{1}$$_{i}$ . These elements, in turn, control dD$_{i}$/dA , Eq. (15). Thus, the progression of W as a function of A can be read off from the plot of D$_{i}$ ( A ), which is given in Fig. 1. It is straightforward to numerically integrate the evolu- tion equations for ( x, y ). To do this we choose the initial values (in W$_{0}$ ) ϵ = 0 . 17, β = 0 . 02, corresponding to the experimental bounds | V$_{e}$$_{3}$ | 2 ≤ 0 . 03 [6] and an assumed CP violation phase cos φ = 1 / 4. Also, ξ = η = 0. In Fig. 2 the results are compared to the approximate solutions obtained earlier (Eqs. (24) and (27)). The agreements are quite good. The evolution of J 2 is shown in Fig. 3. Compared to its vacuum value, it is seen that, except for some enhancement near A = A$_{l}$ , J 2 tends to decrease with increasing A , as one would expect from Eq. (18). In conclusion, in this work we derived the evolution equations for the neutrino parameters as a function of matter density. We found analytic, approximate, but simple solutions of these equations for values centered around the known neutrino parameters. This is possible because of two fortuitous circumstances: 1) the neutrino mass differences are widely separated, enabling one to use the two-flavor resonance approximation; 2) the mixing in 4