A new simulation method for nonuniform polymer solutions between planar surfaces at full chemical equilibrium is described. The technique uses a grid of points in a two-dimensional thermodynamic space, labeled by surface area and surface separations. Free energy differences between these points are determined via Bennett's optimized rates method in the canonical ensemble. Subsequently, loci of constant chemical potential are determined within the grid via simple numerical interpolation. In this way, a series of free energy versus separation curves are determined for a number of different chemical potentials. The method is applied to the case of hard sphere polymers between attractive surfaces, and its veracity is confirmed via comparisons with established alternative simulation techniques, namely, the grand canonical ensemble and isotension ensemble methods. The former method is shown to fail when the degree of polymerization is too large. An interesting interplay between repulsive steric interactions and attractive bridging forces occurs as the surface attraction and bulk monomer density are varied. This behavior is further explored using polymer density functional theory, which is shown to be in good agreement with the simulations. Our results are also discussed in light of recent self-consistent field calculations which correct the original deGennes results for infinitely Iona polymers. In particular, we look at the role of chain ends by investigating the behavior of ring polymers.