Bone is a composite material where tubular osteon structures surrounded by thin interfaces reinforce the tissue at the microscale. These interfaces provide alternative weak paths for crack growth and give rise to potent extrinsic toughening mechanisms. Numerical models can help to gain a deeper understanding of fracture resistance in bone, however, capturing the interface mechanics is key. For this purpose, the phase field method for fracture is appealing due to its capability of capturing complex cracking phenomena. Still, it is far from well-established in the field of biomechanics. In this study, we evaluated a selection of recent open-source implementations of the phase field method to find the best approach for simulating crack growth in bone tissue. We also proposed a new method for correcting the interface fracture toughness in a phase field model. The selected implementations were compared using benchmark tests, including single edge notched tensile specimens with and without interfaces, as well as models representing typical bone geometries. We found the quasi-Newton monolithic solvers to be superior both in terms of computational speed and robustness compared to the evaluated staggered solvers. We also showed that correcting the fracture toughness of the interface is key for simulating crack patterns that are consistent with analytical predictions from linear elastic fracture mechanics.