Abstract More than 50% of proven hydrocarbon reserves are in naturally fractured reservoirs (NFR) and fractures are present in most carbonate reservoirs. The heterogeneity of NFR material properties; permeability and porosity, is of large scale thus presenting a challenge in modeling fluid flow and recovery mechanisms numerically. For a large class of fractured reservoirs, especially for multiphase production mechanisms, the dual-medium approach, using a transfer function (TF) to represent the exchange term between fracture (flowing domain) and matrix (stagnant domain) is a good tool. However, the complexity of modeling the geometry of natural fractures and the advances in finite element method and computing power are calling for newer developments and tools to better model and simulate the unique two-timescale flow-system of these reservoirs. The evaluation of various numerical methods in modeling fluid flow in highly heterogeneous NFR, such as finite element, finite volume, and finite difference, is the main goal of this research. This will be done through writing and benchmarking a simulator for each method and testing its capabilities in capturing the correct recovery mechanisms of NFR while maintaining accuracy in the numerical calculations. In this poster, we first present a literature review followed by an overview of the main equations used in fluid flow in reservoir engineering and their derivation and algebraic approximations in finite difference and finite element methods. After that, we present the results of the first 2D simulator written, using finite element finite volume (FEFV) and our preliminarily results in one and two phase (water-oil incompressible fluid) simple and heterogeneous models. Finally, a future plan of this research is discussed.