Energy systems at district/urban level are getting more complex and diversified from time to time. Different energy carriers are coupled each other to meet various types of energy demands. The conventional way of analyzing energy networks independently does not reflect the true nature of the coupled networks. One of such a promising coupled multi-carrier energy system (MCES) is the combination of district heating and electricity networks. The coupling between these two networks is increasing due to the integration of co- and poly-generation technologies at the distribution networks. Recent literatures tried to address a load flow analysis for lightly coupled networks by formulating case-specific load flow models. This paper presents a more general and flexible tool developed using Matlab (R) which can be used to conduct the load flow analysis of highly coupled electricity and heating networks. An energy hub concept is extended further to formulate a general model in which local generations and detailed network parameters of MCES can be taken into account. Coupled heating and electricity networks are modeled in detail for illustration. The flexibility and generality of the model are then tested by considering case studies with different network topologies (tree and meshed). A comparison is also made with a model used in recent literature. The results show that the proposed model is more accurate. The main contribution of this paper can be summarized by the following five points: (1) Coupling matrices are used to relate network power flow equations of different energy carriers; (2) Hybrid hydraulic head and pipe flow equations are used to develop the hydraulic model which can be applied for both types of tree and meshed heating networks with the possibility of pumping units; (3) A general thermal model that relates steady state temperature drops and mass flow rates, even during change of flow direction, is developed for the heating network; (4) The electricity network is modeled with the possibility of tap changing transformers; (5) The overall system of equations are solved as a single problem using Newton-Raphson iterative method.