Evolution of the mass function (MF) and radial distribution (RD) of the Galactic globular cluster (GC) system is calculated using an advanced and a realistic Fokker–Planck (FP) model that considers dynamical friction, disc/bulge shocks and eccentric cluster orbits. We perform hundreds of FP calculations with different initial cluster conditions, and then search a wide‐parameter space for the best‐fitting initial GC MF and RD that evolves into the observed present‐day Galactic GC MF and RD. By allowing both MF and RD of the initial GC system to vary, which is attempted for the first time in the present Letter, we find that our best‐fitting models have a higher peak mass for a lognormal initial MF and a higher cut‐off mass for a power‐law initial MF than previous estimates, but our initial total masses in GCs, MT,i= 1.5–1.8 × 108 M⊙, are comparable to previous results. Significant findings include that our best‐fitting lognormal MF shifts downward by 0.35 dex during the period of 13 Gyr, and that our power‐law initial MF models well‐fit the observed MF and RD only when the initial MF is truncated at ≳105 M⊙. We also find that our results are insensitive to the initial distribution of orbit eccentricity and inclination, but are rather sensitive to the initial concentration of the clusters and to how the initial tidal radius is defined. If the clusters are assumed to be formed at the apocentre while filling the tidal radius there, MT,i can be as high as 6.9 × 108 M⊙, which amounts to ∼75 per cent of the current mass in the stellar halo.