This study presents a novel approach for analyzing the interaction between two parallel surfaces grafted with polymer brushes in a good solvent. In the proposed approach, molecular dynamics simulations are performed to establish the mean brush height and the standard deviation of the brush height distribution for a given value of the surface separation. The corresponding probability density function (PDF) of the brush height is then determined and a statistical technique is applied to compute the corresponding interaction free energy per unit area of the grafted substrates. Finally, the Derjaguin approximation is employed to determine the corresponding value of the interaction force between the two surfaces. At relatively high surface grafting density as well as under low to moderate compressions of these two parallel plates, the interdigitation effect of the brushes is quite weak and is not considered in the present study. The results obtained for the interaction free energy and interaction force are compared with those derived using the Alexander and de Gennes (AdG) model [1977, “**Adsorption of Chain Molecules With a Polar Head. A Scaling Approach**
,” J. Phys. (Paris), 38, pp. 983–989, 1985, “**Films of Polymer-Solutions**
,” C. R. Acad. Sci., 300, pp. 839–843] and the Milner, Witten, and Cates (MWC) model [1988, “**Theory of the Grafted Polymer Brush**
,” Macromolecules, 21, pp. 2610–2619], respectively. The value of the normalized interaction free energy computed using the present method is higher than that obtained from the AdG and MWC models at larger surface separations. However, the three sets of results are in good agreement particularly at smaller values of the surface separation. In addition, the results obtained by the current method for the interaction force are found to be in better agreement with the experimental data than those obtained using the AdG or MWC models. The enhanced performance of the proposed method is attributed primarily to the use of an adaptive non-Gaussian PDF of the brush height to model the effects of fluctuations in the brush conformation at different distances from the grafting plane.