diff --git a/spatial/r3/vector.go b/spatial/r3/vector.go
index e8e69bd2..f66e5c05 100644
--- a/spatial/r3/vector.go
+++ b/spatial/r3/vector.go
@@ -104,9 +104,9 @@ func Gradient(p, step Vec, field func(Vec) float64) Vec {
 	dy := Vec{Y: step.Y}
 	dz := Vec{Z: step.Z}
 	return Vec{
-		X: field(Add(p, dx)) - field(Sub(p, dx)),
-		Y: field(Add(p, dy)) - field(Sub(p, dy)),
-		Z: field(Add(p, dz)) - field(Sub(p, dz)),
+		X: (field(Add(p, dx)) - field(Sub(p, dx))) / (2 * step.X),
+		Y: (field(Add(p, dy)) - field(Sub(p, dy))) / (2 * step.Y),
+		Z: (field(Add(p, dz)) - field(Sub(p, dz))) / (2 * step.Z),
 	}
 }
 
diff --git a/spatial/r3/vector_test.go b/spatial/r3/vector_test.go
index 7b6ce2e9..938ee9e1 100644
--- a/spatial/r3/vector_test.go
+++ b/spatial/r3/vector_test.go
@@ -306,8 +306,8 @@ func TestDivergence(t *testing.T) {
 
 func TestGradient(t *testing.T) {
 	const (
-		tol = 1e-12
-		h   = 1e-4
+		tol = 1e-6
+		h   = 1e-5
 	)
 	step := Vec{X: h, Y: h, Z: h}
 	rnd := rand.New(rand.NewSource(1))
@@ -316,7 +316,7 @@ func TestGradient(t *testing.T) {
 			p := randomVec(rnd)
 			got := Gradient(p, step, test.field)
 			want := test.gradient(p)
-			if vecApproxEqual(got, want, tol) {
+			if !vecApproxEqual(got, want, tol) {
 				t.Errorf("result out of tolerance. got %v, want %v", got, want)
 			}
 		}