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 naive application of these equations exhibits spurious unphysical properties for even an infinitesimal superfluid gap. We will conclude that it is essential to go {em 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 Poschl-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 Poschl-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 Poschl-Teller models, the most visible consequence of the self-consistent calculation is the change in stability range of the system.