It’s not built-in, if that’s what you’re asking. Might be interesting to do that at some point.
External functions can and have been added to Stan, but it requires non-trivial back end C++ plumbing. Only non-trivial in that it’s a lot of work—it’s all known unknowns (other than finite diff precision).
Finite differences may be old, but I wouldn’t characterize them as good! They’re very slow, requiring on the order of 1 (point and point + epsilon) 2 (point - e, point + e) or 4 (point, point - 2e, point - e, point + e, point + 2e) function evaluations per dimension to use the centered or more stable double versions (as we use, for example, in our testing framework).
They also tend to be low accuracy and unstable near zero. In some cases, you can get enough accuracy in the derivatives to make them usable for HMC. We don’t really know where that boundary is—it certainly depends on the model’s curvature.
Other approaches are to do someting like build an approximation to the black box function, for example using a Gaussian process.