Cohesive fracture, characterized by opening-dependent cohesive tractions in an extended crack tip, is an important phenomenon in metals and plastic materials which generates complex multiscale structures. While phase-field models for brittle fracture have been extensively stud- ied (cf. Ambati et al. [4] and references therein), the first mathematically sound scalar phase- field model was introduced by Conti et al. [1] and was extended only recently in Lammen et al. [2]. The numerical implementation of these emerging models presents significant challenges, particularly in two aspects: solving the resulting strongly coupled system of nonlinear algebraic equations, and managing the large number of degrees of freedom required for accurate simula- tions. To address these challenges, we evaluate the performance of the PUMA Software Toolkit for these problems. This toolkit employs generalized finite-element techniques in a meshfree approach based on the Partition of Unity Method, enabling the use of arbitrary and independent local approximation spaces through enrichment functions. After demonstrating the toolkit’s suitability for these problems, subsequent research work could focus on developing appropriate enrichment strategies — similar to successful applications in peridynamics (Birner et al. [3]) — to substantially reduce the number of degrees of freedom and, consequently, the computa- tional complexity. Additional features of PUMA include parallel computation capabilities and multilevel solvers. We validate the toolkit using several several well-established benchmark tests in two dimensions, covering both brittle and cohesive cases. Our investigation addresses multiple simulation aspects, as the representation of initial cracks, the comparison of staggered and monolithic solving strategies for the coupled system, and the influence of the chosen length scale on the smeared crack model. We benchmark our results against those obtained using the open-source finite-element software FEniCSx.