A variety of building energy analysis and simulation tools are increasingly used to determine peak heating and cooling loads, size thermal plant, anticipate annual energy consumption and analyse thermal comfort. Numerical solution techniques are considered the most flexible for building energy simulation. When applied to the differential equations modelling energy flows in buildings, they give rise to a system of non-linear algebraic (difference) equations. In order to evaluate numerical methods for building energy simulation, the problem has been characterized mathematically and comprehensive test problems (equation sets) with these characteristics have been prepared. The principal attribute of the problem was found to be a stifiess ratio of the order of lo4. Candidate methods have been programmed and their outputs compared, in numerical experiments, with highly accurate (converged) solutions for the test problems. The accepted validation methods, empirical validation, analytical verification and inter-modal comparison were considered inappropriate. The first estimates total and not just numerical error, the second is too confined and the third lacks an absolute standard. The main evaluation parameter used was computational efficiency which is defined as accuracy attained per unit (computational) effort expended. An improved difference equation solver has been proposed and compared with the one used in the European reference model (ESP) and elsewhere. It was found to produce 27% less error than the currently used method. A fundamental method for estimating the pre-conditioning period of a building has been put forward in this part of the work. The trapezoidal rule (TR) is currently used in a number of building energy simulation packages including ESP. A known instability associated with the method is described and an implicit member of the Runge-Kutta family, possessing the necessary strong stability, has been shown, using the test problems, to be more efficient than TR by a factor of 4.27.