Summary: Phylogeny estimation has been the subject of several researches due to its significant importance and numerous applications. The aim of this research is to study the phylogeny estimation from single nucleotide polymorphism (SNP) haplotypes under the maximum parsimony criterion (MPPEP-SNP). Previous exact methods have modeled the mentioned problem as a mixed integer programming (MIP) problem. Since the problem, in general, proved to be NP-hard which causes MIP models to take long runtime for solving large-scale instances, the need for non-exact methods is obvious. In this paper, the authors propose a heuristic algorithm that attempts to find the MPPEP-SNP solution in several stages by solving a specific MIP model in each stage. Created based on network flows formulation, MIP models appearing in each stage are very small; therefore, their exact solutions could be found practically very fast. In order to evaluate the performance of the proposed algorithm, it has been tested on both simulated and real instances and compared with Pars and Flow-RM as two of the best known methods. Our numerical experiments show that the proposed method can compete with the previous methods in terms of accuracy, runtime, and specially the largeness of the solved instances.