In my example above is how I compute the orthogonal components (x,y,z) of the gradient vector. h may even have components itself: hx, hy and hz (e.g. anisotropic grid cells). In this case, it becomes:
void MyFuntion::EvaluateGradient(double* x, double* g)
{
/// hx, hy, hz are small separation values (e.g. grid cell sizes)
g[0] = my_function(x[0] + hx, x[1], x[2]) - my_function(x[0] - hx, x[1], x[2])) / (2 * hx);
g[1] = my_function(x[0], x[1] + hy, x[2]) - my_function(x[0], x[1] - hy, x[2])) / (2 * hy);
g[2] = my_function(x[0], x[1], x[2] + hz) - my_function(x[0], x[1], x[2] - hz)) / (2 * hz);
}