Recent observations of neutron stars with radio, X-rays and gravitational waves have begun to constrain the equation of state for nuclear matter beyond the nuclear saturation density. To one's surprise, there exist approximate universal relations connecting certain bulk properties of neutron stars that are insensitive to the underlying equation of state and having important applications on probing fundamental physics including nuclear and gravitational physics. To date, analytic works on universal relations for realistic neutron stars are lacking, which may lead to a better understanding of the origin of the universality. Here, we focus on the universal relations between the compactness (C), the moment of inertia (I), and the tidal deformability (related to the Love number), and derive analytic, approximate I-Love-C relations. To achieve this, we construct slowlyrotating/tidally-deformed neutron star solutions analytically starting from an extended Tolman VII model that accurately describes non-rotating realistic neutron stars, which allows us to extract the moment of inertia and the tidal deformability on top of the compactness. We solve the field equations analytically by expanding them about the Newtonian limit and keeping up to 6th order in the stellar compactness. Based on these analytic solutions, we can mathematically demonstrate the O(10%) equation-of-state variation in the I-C and Love-C relations and the O(1%) variation in the I-Love relation that have previously been found numerically. Our new analytic relations agree more accurately with numerical results for realistic neutron stars (especially the I-C and Love-C ones) than the analytic relations for constant density stars derived in previous work. Our results provide a mathematical explanation of the amount of universality for realistic neutron stars in the above universal relations.