It does work, but it is even slower than when I define the function as rober(du,u,p,t) and solve it without passing Jacobian. Therefore, I modified my code as follows. The solution is not correct, all the variables are not changing over time.

Note that it automatically uses ForwardDiff.jl. But it’s doing something quite a bit smarter because it’s caching the configs and using SparseDiffTools stuff. If you just want ForwardDiff autodiff’d Jacobians, I’d recommend you just use the default and do nothing.

Hi Chris,
I thought Julia ODE solvers are also using finite difference when Jacobian is not provided by users. Are all the Julia solvers automatically using ForwardDiff.jl and SparseDiffTools? Including CVODE_BDF()?

If I have a very high dimensional and stiff ODE system with no access to the analytical jacobian, the best way to solve it is just to use the default and do not pass a jacobian matrix to the solver. Am I correct?

All of the pure Julia methods default to using ForwardDiff.jl via SparseDiffTools.jl. To turn this off and instead use FiniteDiff.jl, you would set autodiff=false, like Rodas5(autodiff=false). Sundials is an outlier since it uses its native algorithm here. This is documented at Jacobians, Gradients, etc. · DifferentialEquations.jl

I would recommend reading the tutorial on this topic:

Giving a sparsity pattern will automatically specialize the Jacobian calculation on that sparsity pattern via SparseDiffTools.jl, but sometimes analytical solutions to Jacobians can be faster given certain sparsity patterns. That said, using preconditioned Newton-Krylov can also be essential in this domain.