I give here a detailed user guide for the C++ program SHdecay, which has been developed for computing the final spectra of stable particles (protons, photons, LSPs, electrons, neutrinos of the three species and their antiparticles) arising from the decay of a super-heavy X particle. It allows to compute in great detail the complete decay cascade for any given decay mode into particles of the Minimal Supersymmetric Standard Model (MSSM). In particular, it takes into account all interactions of the MSSM during the perturbative cascade (including not only QCD, but also the electroweak and 3rd generation Yukawa interactions), and includes a detailed treatment of the SUSY decay cascade (for a given set of parameters) and of the non-perturbative hadronization process. All these features allow us to ensure energy conservation over the whole cascade up to a numerical accuracy of a few per mille. Yet, this program also allows to restrict the computation to QCD or SUSY-QCD frameworks. I detail the input and output files, describe the role of each part of the program, and include some advice for using it best. Program summary Title of program: SHdecay Catalogue identifier:ADSL Program summary URL: http://cpc.cs.qub.ac.uk/summaries/ADSL Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland Computer and operating system: Program tested on PC running Linux KDE and Suse 8.1 Programming language used: C with STL C++ library and using the standard gnu g++ compiler No. lines in distributed program: 14 955 No. of bytes in distributed program, including test data, etc.: 624 487 Distribution format: tar gzip file Keywords: Super-heavy particles, fragmentation functions, DGLAP equations, supersymmetry, MSSM, UHECR Nature of physical problem: Obtaining the energy spectra of the final stable decay products (protons, photons, electrons, the three species of neutrinos and the LSPs) of a decaying super-heavy X particle, within the framework of the Minimal Supersymmetric Standard Model (MSSM). It can be done numerically by solving the full set of DGLAP equations in the MSSM for the perturbative evolution of the fragmentation functions D p 2 p 1 ( x, Q) of any particle p 1 into any other p 2 ( x is the energy fraction carried by the particle p 2 and Q its virtuality), and by treating properly the different decay cascades of all unstable particles and the final hadronization of quarks and gluons. In order to obtain proper results at very low values of x (up to x∼10 −13), NLO color coherence effects have been included by using the Modified Leading Log Approximation (MLLA). Method of solution: the DGLAP equations are solved by a four order Runge–Kutta method with a fixed step. Typical running time: Around 35 hours for the first run, but the most time consuming sub-programs can be run only once for most applications.