Direct flexoelectricity is a size-dependent phenomenon, very prominent at smaller scales, that connects the strain gradients and the electric field. The very existence of strain gradients enhances noncentrosymmetry and heightens the interaction between piezoelectricity and flexoelectricity, demanding fully coupled higher-order electromechanical formulations. The numerical instability of the existing finite elements used to model flexoelectricity alone is revealed due to their reliance on the stabilization parameter. Thus, two new finite elements () and () are proposed for mixed FEM that are numerically robust without any need of such stabilization parameters. Additionally, the existing finite element [ in (Deng et al. in J Appl Mech 84:081004, 2017)], is implemented from scratch to replicate known results and benchmark the performance of newly proposed finite elements. To verify the robustness of these elements, various benchmark problems for flexoelectricity in dielectric solids, such as a thick cylinder and truncated pyramid are simulated. The great agreement of the numerical results with the existing ones reflects the foundational nature of the proposed elements. Furthermore, the proposed mixed finite elements were used to successfully analyze cantilever beam and truncated pyramid problems where piezoelectric effects were taken into account for the first time. Current results are intrumental in simulating flexoelectricity and piezoelectricity together to highlight their interactions using newly proposed numerically robust finite elements.