Why is it that when I run scipy.integrate.solve_ivp with the Radau method, the number of Jacobian evaluations doesn't equal the number of function evaluations? Specifically, in my code, it only evaluates the Jacobian once or twice, even though the Jacobian changes at every step.
Jacobian evaluation and factoring is an expensive operation. For the correct and efficient working of the ODE integrator an approximate Jacobian is sufficient. The solver algorithms a once computed Jacobian for as long as possible before renewing it. So if the Jacobian change little over the part of the domain of the ODE where the solution segment lives, then one Jacobian evaluation can be sufficient.
Note also that there can be 3,6,9 or 12 or even higher multiples of 3 as number of function evaluations for a single Radau5 step.