Basing on discretizations of Laplace-Beltrami operator and Gaussian curvature over triangular and quadrilateral meshes and their convergence analyses, this paper proposes in this paper a novel approach for constructing geometric partial differential equation (PDE) Bézier surfaces, using several fourth order geometric flows. Both three-sided and four-sided Bézier surface patches are constructed with G1 boundary constraint conditions. Convergence properties of the proposed method are numerically investigated, which justify that the method is effective and mathematically correct.