As far as I can tell, the ODE integrators available in scipy are not constructed with the option to use dask arrays. This imposes an unfortunate limitation when evolving a system of equations on a grid, despite that fact that dask appears an adequate solution. As a practical matter of implementation, how feasible would it be to make this addition to the code base? Perhaps there are a certain set of steps a core maintainer would recommend starting with and someone could take over along those lines if this feature is not given top priority.