Check out the idea of dual numbers, there exist several example implementations in C++. Evaluating the functions using appropriately initialized dual numbers results in the evaluation of one directional derivative. Repeat this for all coordinate directions.
For a nice introduction, see Piponi: AD, C++ and Photogrammetry (http://el.mdu.edu.tw/datacos/09820722022O/paper.pdf)
If L is the effort to evaluate the functions, one directional derivative costs about 3*L, the full jacobian 3n*L. If one were to combine all the directions in one evaluation, this reduces to (1+2n)*L. But by then we have entered the realm of automatic or algorithmic differentiation.
Look for FADBAD/TADIFF for an easy implementation, for really fast codes one uses code transformations as provided by the Tapenade project. ADOL-C is in between, more automatic than Tapenade, faster than FADBAD.