Given a univariate uni-modal density function $f(x)$ (very hard to compute its cumulative distribution function (CDF) $F(x)$, not to mention its inverse CDF $F^{-1}(x)$),
how to find the best Gaussian/normal approximation without drawing samples from the density function $f(x)$ via rejection method and then computing the mean and variance?
There's no need for rejection – you can estimate the mean and variance by numerical quadrature. That only requires some function evaluations; if even that is too expensive, I don't see what you could possibly do, since you'd essentially have a mystery function in a black box that you can't know anything about.