Hot massive stars have strong stellar winds with velocity 1000-3000 km/s that, when shocked, produce hot gas with temperature 107−108 K. This is observed as diffuse X-ray emission in Wolf-Rayet nebulae, around young star clusters, and around the nearest O star, Zeta Oph. Zeta Oph is a runaway star whose strong wind produces a bow shock from the wind-ISM interaction, and diffuse X-ray emission from the shocked wind was detected in Chandra observations. It therefore provides an excellent testbed for models of bow shocks around massive stars. We present the first 3D magnetohydrodynamic simulations of the bow shock around Zeta Oph, calculating the X-ray luminosity and thermal emission maps as a function of energy. Simulations with a range of different assumptions are compared with infrared and re-analysed X-ray observations. In all cases the simulations underpredict the X-ray emission, and the observed morphology of the emission is not easily understood from models. This first numerical study of the bow shock and wind bubble around Zeta Oph uses a relatively simple model of a uniform ISM and ideal-magnetohydrodynamics, and can be used as a basis for comparing results from models incorporating more physical processes. The faint diffuse emission is difficult to characterise from Chandra observations, and the much greater sensitivity to such nebulae promised by ATHENA will dramatically improve our understanding of the hot ISM, stellar winds, and mixing between hot and warm ISM phases.