To accurately analyze the fracture behavior of piezoceramics at small length scales, flexoelectricity must be considered along with piezoelectricity. Due to its dependence on size, flexoelectricity predominates at the micro- and nano-length scales. Additionally, crack tips having the largest strain gradient state cause large flexoelectricity around them. Different approaches are employed in the past to model cracks computationally. However, extended isogeometric analysis (XIGA) is proven to be an accurate and efficient method. C¹ continuity requirements for modeling gradients in flexoelectricity are met by non-uniform rational B-splines (NURBS) basis functions used in XIGA. In this work, XIGA-based computational model is developed and implemented to study the fracture behavior of the piezoelectric-flexoelectric domain. An in-house MATLAB code is developed for the same. Several numerical examples are studied to ensure the efficacy and efficiency of the implemented model, and crack behavior is presented in the form of an electro-mechanical J-integral. The analysis is carried out to investigate how cracks behave for different flexoelectric coefficients under different electrical and mechanical loading combinations. J-integral is also analyzed against crack parameters such as crack orientation and length. It is observed that boundary loads and flexoelectric material constants significantly influence J-integral. Results also show a considerable amount of fracture toughening in the presence of flexoelectricity. The peak value of J-integral is found to be reduced with an increase in the flexoelectric coefficient. A significant reduction in J-integral, as much as 45%, is observed when the flexoelectric constant varied from 0.5 to 2 µCm⁻¹.