Despite the overwhelming focus on the overall average treatment effect in the methodological and statistical literature, in many cases the efficacy of an educational program or intervention might vary based on unit background characteristics. The identification of subgroups for which an educational intervention is particularly effective or, on the other hand, has no effect or is possibly harmful, may have important practical implications, especially in terms of allocation of resources. We propose a five-step approach using propensity score matching and regression trees to identify subgroups with heterogeneous treatment effects in observational studies. Results of two Monte Carlo simulation studies demonstrate that the proposed approach can accurately identify heterogeneous subgroups while maintaining Type I error rate. In a case study with Early Childhood Longitudinal Study-Kindergarten cohort data, we find that the effect of exposure to special education services on fifth-grade mathematics achievement varies based on kindergarten mathematics achievement and student gender.