Bayesian emulation is used in modelling complex simulators and is seen as an efficient and powerful statistical tool. The simulations can be very time-consuming to run, resulting in only being able to run the simulator a limited number of times. From the literature, it is has been suggested that the emulator is best suited for continuous functions; however, it is very common to find physical problems containing discontinuities. These discontinuities' positions may also be unknown and therefore, for this research, information on where the discontinuities will be withheld from the emulator. This thesis focuses on emulating the Heaviside function as one simple function containing step-discontinuities and then progresses to slightly more complex functions with step-discontinuities. Specific goodness-of-fit measure have been designed to highlight and measure the emulator when applied to theses step-discontinuous, such as mean squared error and a simple design of Jaccard index. The numerical calculations of the goodness-of-fit techniques are carried out in the R statistical programming language, with the BACCO package for the emulator. It is found that the emulator is able to model discontinuous functions to some degree but, unless there is certainty about the discontinuities' locations, care should be taken when using Bayesian emulation to model discontinuous functions.