Here, time-dependent solutions for the heat conduction equation are numerically evaluated in 1 D space using a fully implicit algorithm based on the finite difference method, assuming temperature-dependence of thermal conductivity. The method is implemented using the package ‘GEOTHERM’, comprising 13 MATLAB-derived scripts and 3 Excel spreadsheets. In the package, the initial state of the modeled crust, including its thickness, average density, and average heat production rate, can be configured by the user. The exhumation/burial history and metamorphic evolution of the crust are simulated by changing these initial values to fit the vertical displacement rates of the crust imposed by the user. Once the inputs have been made, the variations with depth of temperature, proportion of melt, and shear stress, as well as average values of heat flow at the surface and across the Moho, are calculated and displayed in five separate plots. The code is demonstrated with respect to the Carboniferous evolution of the South Variscan Belt. The best fit to independent petrologic constraints derived from thermobarometry is obtained with an early Carboniferous(342 Ma) slab break-off and a shear strain rate of 10 −13 s −1 between 318 and 305 Ma.