improve computation of epsilon in jacobian_fd

Author: JoshuaLampert

src/semidiscretization/semidiscretization.jl

@@ -240,7 +240,12 @@ function jacobian_fd(semi::AbstractSemidiscretization;
     # use second order finite difference to estimate Jacobian matrix
     for idx in eachindex(u0_ode)
         # determine size of fluctuation
-        epsilon = sqrt(eps(typeof(u0_ode[idx])))
+        # This is the approach used by FiniteDiff.jl to compute the
+        # step size, which assures that the finite difference is accurate
+        # for very small and very large absolute values `u0_ode[idx]`.
+        absstep = sqrt(eps(typeof(u0_ode[idx])))
+        relstep = absstep
+        epsilon = max(relstep * abs(u0_ode[idx]), absstep)
 
         # plus fluctuation
         u_ode[idx] = u0_ode[idx] + epsilon