Invariant differential cross section of $\pi^0$ (left) and $\eta$ (right) versus transverse momentum for pp collisions at $\sqrt{s}=13$ TeV. The data are parametrized with a modified TCM model (see Eq.6.2.) and compared to predictions from PYTHIA8 Monash, PYTHIA8 Ropes, EPOS LHC and predictions from NLO pQCD calculations using recent PDFs and FFs. Ratio plots of the data and model calculations to the modified TCM fit of the data are shown in the lower panels. Statistical error bars are represented by vertical bars, and systematic uncertainties are shown as boxes.