A general scheme is presented to calculate high-temperature series coefficients for ensemble averages of spin operators for spin systems with Hamiltonians containing a large number of model parameters. The scheme, which is based on the moment method, provides the series coefficients as exact functions of the model parameters, e.g., spatial dimensionality, coupling distributions in coordinate and spin space, site-dependent field distributions, and spin quantum number. General expressions for the series coefficients for the auto- and pair-correlation functions are given to sixth order in the case of a classical Hamiltonian with bilinear interactions and a one-component site-dependent magnetic field. The general expressions are used to calculate susceptibility series for the simple cubic anisotropic classical Heisenberg antiferromagnet in a uniform nonordering magnetic field along the easy axis. The series coefficients are polynomials in three variables representing the field, the anisotropy, and the ratio of nearest- and next-nearest-neighbor couplings, respectively. From an analysis of the ordering susceptibility series the phase diagram spanned by the temperature and the field has been calculated for various values of the anisotropy parameter. The calculated phase diagram, which includes a spin-flop phase, an antiferromagnetic phase, and a paramagnetic phase, is in agreement with predictions based on Monte Carlo and renormalization-group calculations.