We present details of a Monte Carlo simulation code which is coupled to a Heat Diffusion Equation (HDE) solver. Through an iterative procedure, which bypasses the differences in electronic and thermal timescales, this coupled code is capable of producing steady-state thermally self-consistent device characteristics. Electronically-generated thermal flux is calculated by monitoring the net rate of phonon emission, which may be resolved both spatially and by phonon type. The thermal solution is extracted through use of a novel analytical thermal resistance matrix technique which avoids calculation of temperatures beyond the electronically important device region while including the large-scale boundary conditions. On application to a GaAs MESFET the expected 'thermal droop' behaviour is obtained in the I-V characteristics and we find a linear relationship between peak lattice temperature and applied source-drain bias. At moderate biases the contribution of intervalley phonons to the thermal power output surpasses that of optical phonons.