A program was written in Python to simulate nuclide reactions and burnup in a thermal fission reactor numerically. The program focused on the depletion calculations and used a simplified neutron flux equation. Nuclide data like cross-sections and fission products were read in from ENDF format files that have undergone pre-processing. To solve the more than 500 simultaneous differential equations that describe the varying isotopic concentrations, short-lived decay isotopes and their decay chains were identified and solved with a modified Bateman solution and then the long-lived isotopes concentrations were solved with matrix exponentiation. The flux was calculated to keep the heat output of the reactions constant. The simulation calculations were validated by comparing the output of decay chains with known analytical solutions. The output of the reactor burnup simulation was compared to that of ORIGEN (The Oak Ridge National Laboratory Isotope Generation And Depletion Code) for a Light Water Reactor at constant load to a burnup of 33GWd/ton. The output of the simulation was relatively similar to that of ORIGEN, but differed in some marked ways, e.g. plutonium breeding, which suggested that the neutron flux calculations and neutron absorption by U238 was not similarly modelled as in ORIGEN. By slightly adjusting the neutron absorption of U238 in the simulation, the correspondence between the simulation and the reference output was improved
