[Roadmap area] Solver improvements

[brosaplanella]

I'd like to continue making pybamm easier to work with for the parameter inference libraries like PyBOP (https://github.com/pybop-team/PyBOP). At the moment I think this includes:

• GPU as well as multithreaded support for running many simulations in parallel (Add support for MLIR-based expression evaluation #3826)
• Integrate sensitivities more into Simulation (generalise simulation, include sensititivies #3834)
• Run multiple simulations from different starting points (IDAKLU solver: add option for multiple initial conditions #3713). Should be possible to reset the starting points during the middle of a solve
• This is generally "making it faster", so profiling would fit in here (see delay xarray.DataArray initialization #3862)
• others?....

Originally posted by @martinjrobins in #3839 (comment)