We present a molecular theory of the phenomena of single-chain collapse and phase separation into a polymer-rich and a polymer-poor phase, which occur in polymer solutions below the Theta temperature. The theory extends the Fourier self-consistent approach of Allegra and Ganazzoli from the study of single-chain properties to that of an ensemble of interacting chains. We derive an expression for the free energy of a "Gaussian cluster" made up of nu chains of length N (nu=1,2,3,...; N much greater than 1). In the limit nu-->infinity this yields a mean-field expression for the solution free energy per chain as a function of the reduced temperature tau(T-Theta)/T, the polymer volume fraction phi and the mean-square radius of gyration of the chains. From this we calculate the chain dimensions in solution and several thermodynamic properties, such as the osmotic pressure and the polymer-solvent phase diagram. We find that the contraction ratio of the chain radius of gyration is a single-valued function of (tau B + K1F phi) root N, where B and K-1 specify the strength of the two- and three-body interactions and F is a polymer-dependent positive constant. We provide numerical evidence for a possible universality of the binodal line for different polymer-solvent systems; the spinodals do not share this characteristic of universality.