Shale gas reservoirs are characterized by gas adsorption and multiple flow mechanisms in shale matrix and complex fracture networks after hydraulic fracturing stimulation. It is of significance to capture the dynamic pressure-dependent properties, including gas PVT properties, gas desorption, flow mechanisms, and stress-dependent fracture permeability. In this paper, this problem is handled by proposing a novel Green element method, which is a good extension of the traditional boundary element method. The domain is discretized into many Cartesian grids, and boundary element method is applied to each element. Discrete fractures are flexibly handled by embedding the fractures into the background grids for the domain, and finite difference method is applied to model flow within the fracture system. The flow between the fracture system and matrix system is coupled to obtain the eventual matrix equation, and Picard successive approximation method is applied to obtain the solution of the nonlinear system. The solution of the novel Green element method is validated by the commercial simulator Eclipse for a fractured tight gas well model. A field case from Longmaxi formation in Southwestern China is used to verify the utility of the model. Moreover, the effects of nonlinear parameters on the performance of the well are analyzed. The results show that the nonlinear factors in shale gas flow model have a great impact on the eventual results of gas production forecasts, especially gas PVT properties and stress-sensitive fracture conductivity.