Recent observations of neutron stars (NS) with radio, X-rays and gravitational waves have begun to constrain the equation of state (EoS) for nuclear matter beyond the nuclear saturation density. To one's surprise, there exist approximate universal relations connecting certain bulk properties of NSs that are insensitive to the underlying EoS and having important applications on probing fundamental physics including nuclear and gravitational physics... To date, analytic works on universal relations for realistic NSs 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), moment of inertia (I), and tidal deformability(related to the Love number), and derive analytic, approximate I-Love-C relations. To achieve this, we construct slowly-rotating/ tidally-deformed NS solutions analytically starting from an extended Tolman VII model that accurately describes non-rotating realistic NSs, allowing 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%) EoS 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 NSs(especially the I-C and Love-C ones) than the analytic relations for constant density stars derived in previous work. Based on these analytic findings, we attribute a possible origin of the universality for the I-C and Love-C relations to the fact that the energy density of realistic NSs can be approximated as a quadratic function, as is the case for the Tolman VII solution. (read more)