We address herein the AT1-version of the phase field method (AT1-PFM) applied to linear elastic heterogeneous materials with continuous variation of Young modulus and fracture toughness \(G_{Ic}\) along crack propagation and crack initiation directions. The method is investigated through both finite element simulations and semi-analytical formulations. In crack propagation scenarios, we show that using a naive \(G_{Ic}\) correction yields results in close agreement with classical Fracture Mechanics predictions, within \(\pm 3\) %, and that energy convergence is achieved as the regularization length \(\ell _0\) decreases, indicating potential \(\Gamma \) -convergence. For crack nucleation, we derive semi-analytical pseudo-homogeneous damage profiles in heterogeneous domains and validate them against finite element results. We find that damage evolves until an ultimate strain beyond which a localized crack forms orthogonally to the material variation, marking a bifurcation from the initial solution. Numerical stability is assessed through the second Gateaux derivative, and the continuation algorithm is used to track the post-instability path. These findings advance the understanding of crack nucleation and propagation in heterogeneous materials and offer a validated modeling framework for future fracture prediction with PFM in Functionally-Graded Materials and human long bones.