We develop the variational and correlated basis functions/parquet-diagram theory of strongly interacting normal and superfluid systems. The first part of this contribution is devoted to highlight the connections between the Euler equations for the Jastrow-Feenberg wave function on the one hand side, and the ring, ladder, and self-energy diagrams of parquet-diagram theory on the other side. We will show that these subsets of Feynman diagrams are contained, in a local approximation, in the variational wave function.In the second part of this work, we derive the fully optimized Fermi-Hypernetted Chain (FHNC-EL) equations for a superfluid system. Close examination of the procedure reveals that the naïve application of these equations exhibits spurious unphysical properties for even an infinitesimal superfluid gap. We will conclude that it is essential to go beyond the usual Jastrow-Feenberg approximation and to include the exact particle-hole propagator to guarantee a physically meaningful theory and the correct stability range.We will then implement this method and apply it to neutron matter and low density Fermi liquids interacting via the Lennard-Jones model interaction and the Pöschl-Teller interaction. While the quantitative changes in the magnitude of the superfluid gap are relatively small, we see a significant difference between applications for neutron matter and the Lennard-Jones and Pöschl-Teller systems. Despite the fact that the gap in neutron matter can be as large as half the Fermi energy, the corrections to the gap are relatively small. In the Lennard-Jones and Pöschl-Teller models, the most visible consequence of the self-consistent calculation is the change in stability range of the system. 7 Discussion 75 Appendix A Cluster expansions for the generating function. Diagonal terms 76 Appendix B Cluster expansions for the generating function. Off-diagonal terms 78 Appendix C Cluster expansions for the energy numerator terms 79 Appendix D Calculation of exchange diagrams 84 1. IntroductionThe repertoire of methods for the quantitative microscopic description of normal quantum many-body systems has condensed, over the past few decades, to a relatively small number of techniques. These can be roughly classified on the one hand side as various shades of large-scale numerical simulation methods and, on the other hand, diagrammatic approaches summing, in some approximation, the parquet class of Feynman diagrams. Numerical simulations are capable of high precision but are computationally demanding and limited to relatively simple Hamiltonians and mostly ground state properties. They also need the physical intuition of the user about the possible state of the system. Semi-analytic, diagrammatic approaches lead to a better understanding of the underlying physical mechanisms, but have limited accuracy due to some necessary approximations. These diagrammatic approaches are versions of Green's function methods [1], Coupled Cluster theory [2], and variational methods [3]. The interconnections between the various me...