In order to investigate pressure performance of multiple fractured horizontal wells (MFHWs) penetrating heterogeneous unconventional reservoir and avoid the high computational cost of numerical simulation, a semi-analytical model for MFHWs combining Green function solution and boundary element method has been obtained, where the reservoir is divided into different homogeneous substructures and coupled at interface boundaries by plane source function in a closed rectangular parallelepiped. Hydraulic fractures are assumed uniform flux and dual porosity model is used for natural fractures system. Then the model is validated by compared with analytical solution of MFHWs in a homogeneous reservoir and trilinear flow model, which shows that this model can achieve high accuracy even with a small interface discretization number, and it can consider the radial flow around each hydraulic fractures. Finally, the pressure responses with heterogeneous parameters of reservoirs are discussed including heterogeneous permeability, non-uniform block-length and fracture half-length distribution as well as dual porosity parameters like elastic storage ratio and crossflow ratio.