The covariant variational approach to Yang-Mills theory is further developed. After discussing the foundations of the method both at zero and finite temperature, we briefly recall the effective action for the Polyakov loop and the critical properties of the deconfinement phase transition within this approach. The thermodynamics of pure Yang-Mills theory are studied in detail and the resulting equation of state is compared to lattice data. While there is good agreement in the deconfined (high-temperature) phase, a small but non-zero pressure is predicted in the confined phase at low temperatures, in contrast to physical expectations. We propose possible improvements to address this issue. Finally, we discuss the combination of the variational approach with Dyson-Schwinger techniques and argue that the method can be used as a tool to determine the optimal vertices for a truncated set of Dyson-Schwinger equations. We briefly lay out how this technique could be applied to Yang-Mills theory at zero temperature.