The Cryogenic Underground Observatory for Rare Events (CUORE) is the most sensitive experiment searching for neutrinoless double-beta decay (0νββ) in 130Te. CUORE uses a cryogenic array of 988 TeO2 calorimeters operated at ∼10 mK with a total mass of 741 kg. To further increase the sensitivity, the detector response must be well understood. Here, we present a non-linear thermal model for the CUORE experiment on a detector-by-detector basis. We have examined both equilibrium and dynamic electro-thermal models of detectors by numerically fitting non-linear differential equations to the detector data of a subset of CUORE channels which are well characterized and representative of all channels. We demonstrate that the hot-electron effect and electric-field dependence of resistance in NTD-Ge thermistors alone are inadequate to describe our detectors' energy-dependent pulse shapes. We introduce an empirical second-order correction factor in the exponential temperature dependence of the thermistor, which produces excellent agreement with energy-dependent pulse shape data up to 6 MeV. We also present a noise analysis using the fitted thermal parameters and show that the intrinsic thermal noise is negligible compared to the observed noise for our detectors.