The Galla system:

a simulation study of the case of multiple sons.

Olinick describes a deterministic model first, due to Hoffmann[1965] (http://ceadserv1.nku.edu/longa/classes/mat375/Markov/Hoffmann-1965.pdf)

Let’s use these two distributions for the simulation:
a poisson (to decide how many children), and
a normal (to choose when they were born).

λ=3 (* this is the mean for the poisson: *)
DiscretePlot[PDF[PoissonDistribution[λ],k]//Evaluate,{k,0,10},PlotRange->All,PlotMarkers->Automatic]
σ=5; (* this will be the standard deviation of the son birthday distribution *)
Plot[PDF[NormalDistribution[40,σ],x],{x,20,60}]

We can use them to choose the number of sons that a man has, and then their birthday relative to his age (with mean 40):

(* how many sons does a father have? *)
howManySons[]:=RandomVariate[PoissonDistribution[3]]
howManySons[]
haveASon[number_]:=
Table[RandomVariate[NormalDistribution[40,σ]],{i,1,number}]
haveASon[3]
haveASon[0]

Now we create a “generations” function, that starts with a single father’s age at entry, and pushes his children forward for the stated number of furtherGenerations:

generations[parentAgeAtEntry_, furtherGenerations_]:=
If[furtherGenerations==0, (* We've come to the end of the line: *)
parentAgeAtEntry,
Module[
{sons=howManySons[]},
If [sons ==0,  (* Simply report the would-be parent's age: *)
parentAgeAtEntry,
(* otherwise, recurse, dropping down one further generation: *)
Module[
{
sonsPs=haveASon[sons], (* this determines when they were born, relative to the father *)
sonsAges
},
sonsAges=parentAgeAtEntry+ 40 - sonsPs; (* formula from Olinick *)
Flatten[Table[generations[sonsAges[[i]],furtherGenerations-1],{i,Length[sonsAges]}]]
]
]
]
]

And now we  actually compute some generations, and have a look at the age distribution of the sons at their time of entry into the grade system:

agesAtEntry=generations[20,8];
StandardDeviation[agesAtEntry]
Length[agesAtEntry]
Mean[agesAtEntry]
Histogram[agesAtEntry]