I would like to use ForwardDiff.jacobian to provide Jacobian matrix to ODE solvers. I wonder if it can be faster than using numerical Jacobian. I think the rhs of differential equations should in form of f(du, u,p,t) and automatic jacobian should looks like ForwardDfiff.jacobian(f,u). But I am not sure how to make these two format working together. My code is attached.