The deep coalbed methane (CBM) development is controlled by high temperature, high pressure, and strong in-situ stress, and the coupling effects among the temperature field, stress field, and seepage field are significant. After hydraulic fracturing, reservoir heterogeneity around the fractured well is strong. These factors make it difficult to accurately characterize the complex non-steady-state flow behavior, and thereby increase the uncertainty in predicting the well productivity. Here, a fully coupled mathematical model of non-steady gas-water two-phase flow in heterogeneous formation under the coupling of temperature field, stress field and seepage fields is established for deep coalbed methane development by fractured vertical wells. Moreover, a Gaussian distribution function is employed to describe the spatial heterogeneity of formation permeability around the fractured vertical well. The coupled multi-physical model is numerically solved by using the finite element software COMSOL, and its validity is also verified through numerical simulation history matching of field production data of a deep coalbed methane well at Yichuan area in Ordos Basin. Influencing factor analysis on productivity indicates that the multi-field coupling effect has a significant impact on the prediction of early daily production and late cumulative production, and ignoring any physical field will lead to large computation deviation; for the multi-physical field coupling model, ignoring the temperature field can lead to an approximately 8% underestimation of EUR, while neglecting the stress field may result in a 5% overestimation of EUR. It is recommended to prevent hydraulic fracturing energy from being excessively consumed in the near-wellbore zone and maximize the spread range of fracturing energy so as to reduce formation heterogeneity after fracturing, which can effectively expand the propagation range and velocity of pressure drop and adsorbed gas desorption area. This is conducive to the improvement of coalbed methane production and estimated ultimate recovery. All in all, the two-phase seepage flow model of deep coalbed methane wells considering the coupling of temperature field, stress field and seepage field and permeability heterogeneity of fractured reservoirs can accurately predict dynamic production performance after fracturing. It provides theoretical models and technical supports for production allocation, stimulation optimization and EUR evaluation of deep CBM wells.