The Monte Carlo method based on two‐dimensional entropic sampling within the Wang–Landau (WL) algorithm is applied to simulation of a continuous model of a polyelectrolyte between membrane surfaces. Membranes are presented by parallel plane surfaces holding either fixed or mobile dipoles (representing lipid headgroups). A strongly charged polyion accompanied by neutralizing counterions is placed between the membranes. Periodic boundary conditions are imposed along X‐ and Y‐axes. The volume of the main cell is varied during the simulation by shifting one of the surfaces along Z‐axis. Within two‐dimensional WL sampling algorithm we obtain joint density of states as a function of energy and volume in a single run. In order to increase efficiency of our calculations we introduce a number of modifications to the original WL‐approach. Various properties of the system over wide temperature and volume or pressure ranges, i.e., conformational energy, heat capacity, and free energy, are obtained from the two‐dimensional density of states by simple integration. The osmotic pressure is calculated as a derivative of Helmholtz free energy. Alternatively, properties of the system, including average volume, can be obtained under condition of NPT ensemble. It is shown that the both approaches produce coinciding Posm(V) isotherms. In all considered cases we observe repulsive effective interaction between the membrane surfaces, and repulsion is stronger for the surfaces with fixed dipoles.