The (perhaps naive) approach I've taken in the past is to search the space of ODE systems using genetic programming with a loss function obtained by integrating my candidate ODE system and measuring how closely the output matches the data I want to perform regression on.

What they're doing here is not that. I think if I understand it correctly they're using the neural network to generate data which when symbolically regressed with PySR yield the RHS of each ODE in the system.

What's not immediately obvious is what the benefit of introducing the neural network is--does it make it faster than the "direct" naive approach?