I'm using the solve_ivp module in scipy. The LSODA option offers a Python wrapper to a Fortran integration package. The issue is that the Python wrapper integrates using LSODA timestep-by-timestep. So, there can be a lot of Python overhead to the integration if there are many timesteps. I would like to avoid this overhead by integrating it all at once. Any thoughts on how to do this? Looking into the scipy source code, I see that the Python wrapper feeds LSODA a flag itask = 5 which requests integration by one timestep, and I see that there is an alternate flag in the LSODA documentation for itask = 1 indicating integration all the way up to a final time (as many timesteps as are required). The latter seems better for my purposes because it avoids significant Python overhead. The problem is, when I just change the value of itask in the scipy source code, everything breaks. I think the Python wrapper is built for one timestep at a time... Any suggestions? If going through a modified version of the Python wrapper is too difficult, I'd also be happy to call the lsoda routine directly myself. Just don't really know how to do that--I've only ever gone through solve_ivp. Continue reading...