If your function is badly behaved (e.g. noisy, very oscillatory), no method will perform properly (differentiation is numerically very unstable). That being said, for "nice functions", I have good experience with polynomial (Richardson) extrapolation methods. This paper and this paper give hints on how you might write your own implementation. I will note that this is the method implemented in the NAG numerical libraries (with of course a few wrinkles of their own).
There are two possible alternatives if for some reason you don't want to use Richardsonian methods. One is to use Cauchy's differentiation formula:
$$f^\prime(x)=\frac1{2\pi i}\oint_\gamma \frac{f(t)}{(t-x)^2}\mathrm dt$$
where it is up to you to choose a suitable counterclockwise contour $\gamma$ (a circle is customary); the other is to use the "Lanczos derivative":
$$f^\prime(x)=\lim_{h\to 0}\frac{3}{2h^3}\int_{-h}^h t\;f(x+t)\mathrm dt$$
where you either will have to experiment with an appropriate step size $h$, or use some extrapolative procedure.
You will have to experiment with your computing environment to choose.
