h i g h l i g h t s• A symmetry conserving approximation based on Dyson-Schwinger equations is proposed.• The approximation keeps all the Ward Identities in a theory with continuous symmetry. • The approximation is more accurate than other mean-field type approximations.• The penetration depth of Ginzburg-Landau model shows a downward cusp near T c . a b s t r a c tCondensed matter systems undergoing second order transition away from the critical fluctuation region are usually described sufficiently well by the mean field approximation. The critical fluctuation region, determined by the Ginzburg criterion, |T /T c − 1| ≪ Gi, is narrow even in high T c superconductors and has universal features well captured by the renormalization group method. However recent experiments on magnetization, conductivity and Nernst effect suggest that fluctuations effects are large in a wider region both above and below T c . In particular some ''pseudogap'' phenomena and strong renormalization of the mean field critical temperature T mf can be interpreted as strong fluctuations effects that are nonperturbative (cannot be accounted for by ''gaussian fluctuations''). The physics in a broader region therefore requires more accurate approach. Self consistent methods are generally ''non-conserving'' in the sense that the Ward identities are not obeyed. This is especially detrimental in the symmetry * 229 broken phase where, for example, Goldstone bosons become massive. Covariant gaussian approximation remedies these problems. The Green's functions obey all the Ward identities and describe the fluctuations much better. The results for the order parameter correlator and magnetic penetration depth of the Ginzburg-Landau model of superconductivity are compared with both Monte Carlo simulations and experiments in high T c cuprates.