There is indeed an analytic approach.

If a generation contains $n people, then what the script does is adding $n random variables. Since all 6 possible values are equally likely, the resulting sum is a Binomial distribution.

For big $n (let's say $n >= 20) the Binomial distribution can be very well approximated by a normal distribution.

Normally distributed random variables can be easily and efficiently implemented with the Box–Muller transform, which turns two uniformly distributed random numbers (what rand returns) into two normally distributed random numbers.

If this were a project of mine (and I needed an efficient approach), I'd follow the approach of individual random numbers per population item as long as the population is less than 20, and for higher numbers I'd use the normal distribution approximation.

If there's interest, I can also show how the mean and variance of the random variable is calculated, but I'm tired right now and I'd do it tomorrow :-)

Update: mean and variance isn't as complicated as I thought.

The mean is just (0 + 1 + 2 + 3 + 4 + 5) / 6 = 2.5, and the variance is 1/6 * ((0-2.5)**2 + (1 - 2.5)**2 + (2-2.5)**2 + (3-2.5)**2 + (4-2.5)**2 + (5-2.5)**2) = 35/12 = 2.91666666666667. With these parameters you can just use the formulas in the various Wikipedia entries.