In general, oil reservoirs may consist of composite sedimentary structures composed of materials such as sand, clay, or limestone, which exhibit varying lithology due to sedimentary processes. A comprehensive knowledge of this lithology is essential for accurately assessing their hydrocarbon storage and production capacity. Additionally, this information is indispensable for the implementation of various recovery techniques such as waterflooding, gas injection, surfactant injection, polymer injection, alkaline water solution injection, and others. Highly viscous oil can exhibit non-Newtonian behavior during water injection in certain cases. The use of surfactants, alkaline, or polymer solutions in enhanced oil recovery also introduces non-Newtonian behavior. Recovery methods face challenges when non-Newtonian phases, gravity, and reservoir heterogeneity are combined. Against this backdrop, this work presents a mathematical model for immiscible two-phase flow in an axially composite reservoir with a periodic-layered structure. The model considers a non-Newtonian plastic Bingham-type phase extension of the Buckley-Leverett model. To address the porous medium’s heterogeneity with discontinuous flux functions, the numerical solutions were obtained using the Lax-Friedrichs and Lagrangian-Eulerian schemes. The numerical solutions were compared to analytical solutions obtained using an extended version of Oleinik’s geometric construction for discontinuous flux functions. The outcomes display shock and rarefaction waves, as well as a fixed shock due to the porous medium heterogeneity. The numerical results closely correspond with the analytical solutions, seen particularly in the greater accuracy of the Lagrangian-Eulerian method compared to the Lax-Friedrichs method.