Consistent nonlinear extension of the superparamagnetic response theory is constructed on the basis of the rotary diffusion equation governing the precession and relaxation of the magnetic moment in a single-domain particle in the presence of thermal fluctuations. Along with the linear one, we determine the third- and fifth-order susceptibilities to an arbitrarily directed field by a numerical procedure that takes into account the full set of the magnetic moment relaxation modes, i.e., fast intrawell and slow interwell ones. To get to real samples, the averagings over the particle volumes and anisotropy axes distribution are performed. The results were quantitatively tested on the measurements of the cubic susceptibility of nanocrystalline CuCo alloys, where more simple models miss the goal. Also, for the mentioned situations, using numerical solutions as an quality check, we construct compact and accurate semi-heuristic formulas very feasible for experiment interpretation.Another analytical tool that we derive and that is very handy for theory and experiment is the high-barrier approximation. It works when the ratio S of the magnetic anisotropy energy to temperature, is high. Then the magnetic relaxation is rather slow, and one may say that the particles are becoming quasi-permanent magnets. When treated numerically, this case is the most resource-consuming one. Expanding the micromagnetic equation with respect to 1/S, we get the solutions in the form of asymptotic series, of which the set of dynamic susceptibilities is constructed. The emerging expressions appear to be both very simple and very accurate in the S > 1 range. Partial financial support from Award PE-009-0 from CRDF and Project 97-31311 from INTAS is gratefully acknowledged.