Recently, metasurfaces have emerged as a powerful and versatile paradigm for making many optical devices especially optical lenses. However, the approach usually used to engineer metasurface devices assumes that neighboring elements are identical, by extracting the phase information from simulations with periodic boundaries, or that near-field coupling between particles is negligible, by extracting the phase from single particle simulations. This is not the case in general and the approach thus prevents the optimization of devices that operate away from their optimum. We propose a versatile numerical method that we call local phase method (LPM) to obtain the phase of each element within the metasurface while accounting for near-field coupling with different neighbors. The proposed local phase method paves the way to highly efficient metasurface devices including, but not limited to, deflectors, high numerical aperture concentrators, lenses, cloaks, and modulators.