Hadron diffusion equations with energy-dependent interaction mean free paths and inelasticities are solved using the Mellin transform. Instead of using operators on the finite difference terms, the Mellin transformed equations are expanded in a Taylor series to a first-order partial differential equation in atmospheric depth t and in transform parameter s. These equations are then solved by the method of characteristics. The hadron fluxes (nucleon and pion) in real space are evaluated by the method of residues. For the case of a regularized power law primary spectrum, these hadron fluxes are given by simple residues and one, never before mentioned, essential residue. Comparisons of our solutions are made with the nucleon flux measured at sea level and with the hadron fluxes measured at t = 840 g cm−2 and at sea level. The integral hadron flux is also compared with experimental data from the Mt Fuji Collaboration (t = 650 g cm−2). The agreement between all of them is very good in general, better than 90%.