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][1] and [this paper][2] give hints on how you might write your own implementation. I will note that this is the method implemented in the [NAG numerical libraries][3] (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"][4]:

$$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.

  [1]: http://dx.doi.org/10.1007/BF02166671
  [2]: http://dx.doi.org/10.1007/BF01933664
  [3]: http://nag.co.uk/numeric/numerical_libraries.asp
  [4]: http://www.jstor.org/pss/2589707