* We describe the implementation of XFEM in Matlab with emphasis on the design of the code to enable fast and efficient computational times of fracture problems involving many cracks and arbitrary crack intersections. The efficiency and robustness of the implementations are tested via the solutions to several benchmark problems involving multi-crack growth with coalescence. * We carry out verification studies of the implementation of the minimum energy criterion and give comparisons with the maximum tension criterion. The comparisons reveal that both criteria converge to virtually the same fracture solutions albeit from opposite directions. In other words, it is found that the converged fracture path tends to lie in between those obtained by each criterion on coarser meshes. * From the convergence trends observed in the numerical solutions by different crack growth criteria, we propose a numerical improvement to the crack growth direction criterion which is to assume the average direction of the directions obtained by the maximum tension and the minimum energy criteria. The fracture paths obtained by the modified criterion show significant improvements in accuracy and convergence rates of the fracture paths, as observed from a multitude of test cases. * Finally, we make available the open-source Matlab code including documentation, benchmark and example cases as supplementary material for the verification of our results and with the hope that the code can be useful in general.