The accurate modeling and analysis of integrated energy systems (IES) can precisely obtain the system state variation trend and operation characteristics of IES, and provide the foundation and preliminary for the planning, control, and schedule of IES. However, the commonly deployed Newton–Raphson method suffers from the difficulty of initial state setting, and the differences in the operating characteristics of heterogeneous energy systems bring challenges to the modeling and analysis of IES. In this article, a hybrid physics and data-driven approach is proposed along with a thermal unified matrix model for the steady-state and transient analysis of electricity–heat IES (EH-IES). Both integrated and decomposed solution algorithms are developed to obtain the final system states. Numerical results of several EH-IESs demonstrate the effectiveness of the proposed approach, and the convergence and computational efficiency of the state analysis of EH-IES are remarkably improved.