What is the best way to calculate the numerical gradient of a matrix, representing multiple points in 3D space?

If that’s the case then this is probably the simplest answer (modified from this)

function central_difference(v::AbstractMatrix)
    dv = diff(v, dims=1) / 2
    a = [dv[[1], :]; dv]
    a .+= [dv; dv[[end], :]]
    a
end