Add Automatic Differentiation for ODE solvers

[jespfra]

In Julia, you can use automatic Differentiation for the ODE solvers. However, this is so far not possible in this code. The reason is that the code has matrix products with the state. There is a work around that, see (https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#I-get-Dual-number-errors-when-I-solve-my-ODE-with-Rosenbrock-or-SDIRK-methods). However, this requires changing a fair bit of code. Also, you cannot do that with sub-arrays i.e., x[idx]. Hence, it is not prioritized at the moment. It it however stated that the performance decrease for using finite difference is small for BDF solvers which is used in the code.