  The self-consistent Gaussian approximation (SCGA) for classical spin systems
described by a completely anisotropic D-component vector model is proposed,
which takes into account fluctuations of the molecular field and thus is a next
step beyond the molecular field approximation. The SCGA is sensitive to the
lattice dimension and structure and to the form of spin interactions and yields
rather accurate values of the field-dependent magnetization m(H,T) and other
thermodynamic functions in the whole plane (H,T) excluding the vicinity of the
critical point (0,T_c), where the SCGA breaks down, showing a first-order phase
transition. The values of T_c themselves can be determined in the SCGA with an
accuracy better than 1% for actual 3-dimensional structures. At low and high
temperatures the SCGA recovers the leading terms of the spin-wave theory, the
low- and high-temperature series expansions, respectively. The accuracy of the
SCGA increases with the increase of the spin dimension D, and in the limit D
\to \infty the exact solution for the spherical model is recovered.
