I already figured out the issue. Had to replace all of the np.empty inside the Heun function with np.full. The initial conditions were also way too big for the model I'm studying, so had to change that as well. I now get a proper, realistic looking orbit, like this: Anyways, thank you all who tried to help me.