A numerical formulation coupling finite and boundary element methods is developed to analyze the mechanical deformation and electric polarization of flexoelectric solids experiencing geometrically nonlinear deformation. The proposed method considers the electrical interactions among flexoelectric solids, electric charges, and their surrounding medium. First, a higher-order gradient theory is proposed based on the skew-symmetric couple-stress model to analyze the geometrically nonlinear deformation of flexoelectric solids. This theory includes a total Lagrangian weak form that satisfies linear momentum conservation, angular momentum conservation, and Gauss’s law. Based on the proposed theory, a finite element is developed using basis functions that satisfy C1 continuity. Second, a coupled formulation is developed to consider the electrical interactions among solids, electric charges, and their surrounding medium. In this formulation, conventional boundary elements are adopted to account for the electrostatic surroundings. Besides, electric boundary conditions are naturally imposed on solid boundaries according to the electrical interactions between solids and their electrostatic surroundings. Finally, the proposed method is validated via the comparisons of the numerical results with closed-form solutions.