I'm trying to implement a Rosenbrock method, following the procedure laid out in the book Numerical Recipes in C. My problem is that my system is too complicated to determine an analytical Jacobian. I solved the system in Matlab, which did not require inputting a Jacobian as it automatically calculates it using finite differences. I have been trying to figure out how to do it myself. My main problem is the step size $\Delta y$ to be used for calculating the Jacobian. Surely it must be consistent with the $\Delta y$ being used in the ODE solver? Then do I have to iterate between the ODE solver and the Jacobian calculation for determining an appropriate $\Delta y$?