Strongly correlated electron systems are generally described by tight-binding lattice Hamiltonians with strong local (onsite) interactions, the most popular being the Hubbard model. Although the half-filled Hubbard model can be simulated by Monte Carlo (MC), physically more interesting cases beyond half-filling are plagued by the sign problem. One therefore should resort to other methods. It was demonstrated recently that a systematic truncation of the set of Dyson-Schwinger equations for correlators of the Hubbard, supplemented by a "covariant" calculation of correlators leads to a convergent series of approximants. The covariance preserves all the Ward identities among correlators describing various condensed matter probes. While first-order (classical), second-order (Hartree-Fock or Gaussian), and third-order (Cubic) covariant approximation were worked out, the fourth-order (quartic) seems too complicated to be effectively calculable in fermionic systems. It turns out that the complexity of the quartic calculation in local interaction models,is manageable computationally. The quartic (Bethe-Salpeter-type) approximation is especially important in 1D and 2D models in which the symmetry-broken state does not exists (the Mermin-Wagner theorem), although strong fluctuations dominate the physics at strong coupling. Unlike the lower-order approximations, it respects the Mermin-Wagner theorem. The scheme is tested and exemplified on the single-band 1D and 2D Hubbard model.