Use auto-differentiation. Automatic differentiation is faster than other forms of differentiation and gets errors at machine epsilon. However, it's much harder to implement, so you might need to use a package for it. That said, since it's so useful there are plenty of packages, one I use very often is ForwardDiff.jl.
It's a shame that less academics don't use automatic differentiation because it really is what should be the workhorse, and it's just not taught in graduate school for some reason (though those who go into machine learning know it as "backwards propagation").
One can autodifferentiate functions from libraries in Julia since the libraries themselves are written in Julia and are "duck-typed", and so ForwardDiff.jl can auto-differentiate library functions (including ones which use loops, conditionals, etc.) due to its type system. The only problem it runs into is when there are C/Fortran libraries called, so really you just have to avoid library functions which internally call BLAS (so no matrix multiplications and you're good).
ForwardDiff will now use generic linear algebra fallback function to compute the derivatives, so those will be fine as well. The only Base library functions I can think of where automatic differentiation will fail now (in Julia) is when trying to differentiate a function which calculates an FFT, since that uses FFTW and doesn't have a Julia fallback (for now).