Flow in formations having both low-porosity and low-permeability prevalently occurs in fracture networks. Fracture connectivity and permeability are the main features governing flow in these media, and affect the effectiveness of any subsurface operation (e.g. enhanced oil recovery, enhanced geothermal systems). Moreover, the rheology of the fluid is often non-Newtonian due to their complex microstructure, resulting in a non-linear constitutive law at the continuum scale. Modeling flow of complex fluids in a fracture is challenging, when both medium heterogeneity and a realistic fluid rheology, typically shear-thinning, are considered. Full 3-D simulations are computationally expensive and time-consuming; alternatively, a lubrication-based approach can be adopted to implement an efficient 2-D flow solver able to produce vast statistics in reasonable time. We propose a numerical code that solves the 2-D generalized non-linear Reynolds equation for Ellis fluid flow in a single variable aperture fracture, adopting the finite volume method. The inexact Newton-Krylov method solves the associated non-linear system of equations, while a preconditioned conjugate gradient algorithm is used to solve the linear system at each Newton iteration. A parameter continuation strategy is also introduced to handle strongly non-linear cases. A synthetic fracture aperture field w(x) is generated and the flow problem is solved. The numerical scheme is particularly efficient and robust for most input parameters of practical interest. Numerical parameters of the Inexact Newton-Krylov method have been optimized to simulate flow on large mesh (e.g. 10¹⁰x10¹⁰), considering a range of fracture closure (σ_w/‹w›) from 0 to 1, a shear-thinning index (n) of the Ellis rheologic model varying from 0.1 (strongly shear-thinning behavior) to 1 (Newtonian case). Pressure gradients (∇P) typical of subsurface flow, either natural or forced, have been considered to study natural phenomena and industrial applications.