  The nonlinear ac stationary response of uniaxial paramagnets and
superparamagnets - nanoscale solids or clasters with spin number S ~ 10^0 -
10^4 - in superimposed uniform ac and dc bias magnetic fields of arbitrary
strength, each applied along the easy axis of magnetization, is determined by
solving the evolution equation for the reduced density matrix represented as a
finite set of three-term differential-recurrence relations for its diagonal
matrix elements. The various harmonic components of the magnetization, dynamic
magnetic hysteresis loops, etc. are then evaluated via matrix continued
fractions indicating a pronounced dependence of the nonlinear response on S
arising from the quantum spin dynamics. In the linear response approximation,
the results concur with existing solutions.
