In this thesis, a predictive non-intrusive reduced order model of a high-temperature gas-cooled reactor has been created. We first set out to create a representative multi-physics model coupling the neutronics and thermodynamics in the HTGR. To this end we first made a representative model of neutronics and thermodynamics separately in the HTGR. These representative models were used to train and test our reduced order models (ROMs). It has been found that the representative model converges spatially and temporally to the expected values of analytical and pseudo analytical solutions for the neutronics, the temperature, and the coupled model. However, the convergence in the spatial domain was found to be slower than the mathematical expectation from the implemented finite volume method, converging at a rate of Δx1.5, instead of Δx2. The temporal convergence rate was Δt as is expected from the implemented first order backward differentiation formula. In literature review, we have found that Operator Inference provides a way to construct ROMs that are both predictive and non-intrusive. We used Operator Inference to create three ROMs. A ROM of the temperature model, a ROM of the neutronics model, and lastly a ROM of the coupled model. The temperature ROM created with 12 modi and 1500 snapshots was found to be accurate to the representative model for any sort of power input, achieving a root mean square error below 10−5 compared to the representative model, while predicting 3000 seconds longer than the training at Δt = 1 s. for 1000 random power inputs not included in training. The neutronics ROM created with 7 modi and 3960 snapshots retrieved accurate results, with a root mean square error around 10−4 compared to the representative model if we used a homogeneous temperature in the reactor. If we defined more temperature zones in the reactor we used a different interpolation technique. For this we found the results to be of lesser quality with an average root mean square error closer to 10−2 for a temperature profile far from a steady state temperature profile, while the results were of greater quality with an average root mean square error of 10−5 for temperature profiles close to a steady state temperature profile. Finally, we developed two methods of creating a reduced order model of the coupled model. First, we created a sequential model, where the outputs of the individual temperature ROM and neutronics ROM fed into each other. Second, we created a merged model, in which both the temperature and neutronics were determined by our reduced order model in one single model. A loss of forced cooling was simulated by changing the heat transfer coefficient at the coolant channel boundary. We trained the ROMs for 50 seconds with Δt = 10−3 s, with 5 different heat transfer coefficient equidistant between 0 and 0.03 Wcm−2K−1. It has been found that the merged model was robust and provided accurate results, it could predict to 200 seconds estimating for a heat transfer coefficient of 0.00375 Wcm−2K−1. with an average root mean square error in time for both the neutron flux profile and the temperature profile of less than 10−2. The sequential model was unstable, with root mean square errors orders of magnitude higher than the merged model.