The Kugel-Khomskii (KK) Hamiltonian describes spin and orbital superexchange interactions between d(1) ions in an ideal cubic perovskite structure, in which the three t(2g) orbitals are degenerate in energy and electron hopping is constrained by cubic site symmetry. In this paper we implement a variational approach to mean-field theory in which each site i has its own nxn single-site density matrix rho(i), where n, the number of allowed single-particle states, is 6 (3 orbital times 2 spin states). The variational free energy from this 35 parameter density matrix is shown to exhibit the unusual symmetries noted previously, which lead to a wave-vector-dependent susceptibility for spins in alpha orbitals which is dispersionless in the q(alpha) direction. Thus, for the cubic KK model itself, mean-field theory does not provide wavevector "selection," in agreement with rigorous symmetry arguments. We consider the effect of including various perturbations. When spin-orbit interactions are introduced, the susceptibility has dispersion in all directions in q space, but the resulting antiferromagnetic mean-field state is degenerate with respect to global rotation of the staggered spin, implying that the spin-wave spectrum is gapless. This possibly surprising conclusion is also consistent with rigorous symmetry arguments. When next-nearest-neighbor hopping is included, staggered moments of all orbitals appear, but the sum of these moments is zero, yielding an exotic state with long-range order without long-range spin order. The effect of a Hund's rule coupling of sufficient strength is to produce a state with orbital order.