  The variational method and the Hamiltonian formalism of QCD are used to
derive relativistic, momentum space integral equations for a quark-antiquark
system with an arbitrary number of gluons present. As a first step, the
resulting infinite chain of coupled equations is solved in the nonrelativistic
limit by an approximate decoupling method. Comparison with experiment allows us
to fix the quark mass and coupling constant, allowing for the calculation of
the spectra of massive systems such as charmonium and bottomonium. Studying the
results with and without the nonAbelian terms, we find that the presence of the
nonAbelian factors yields better agreement with the experimental spectra.
