We present results from two different methods for the calculation of hadron spectra in QCD and supersymmetric QCD with large primary energies $\sqrt{s}$ up to ${10}^{16}\mathrm{GeV}.$ The two methods considered are a Monte Carlo (MC) simulation and the evolution of fragmentation functions described by the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations. We find that the pion, nucleon, and all-hadron spectra calculated with the two methods agree well. The MC simulation is performed with new hadronization functions (in comparison with our previous work), motivated by low energy $(\sqrt{s}<{M}_{Z})$ data and the DGLAP equation. The hadron spectra calculated with both sets of hadronization functions agree well, which indicates that our method for calculating the hadronization function works successfully. The small difference in the calculated hadron spectra characterizes the uncertainties of this method. We calculate also the spectra of photons, neutrinos, and nucleons and compare them with other published results. The agreement is good for all x from $\ensuremath{\sim}{10}^{\ensuremath{-}5}$ up to $x<~0.3.$ The consistency of the spectra calculated by different methods allows us to consider the spectral shape as a signature of models with decays or annihilations of superheavy particles, such as topological defects or superheavy dark matter. The ultrahigh energy cosmic ray spectra from these sources are calculated.