diff --git a/spatial/r3/mat_test.go b/spatial/r3/mat_test.go
index c6dd58c9..22a1a5ee 100644
--- a/spatial/r3/mat_test.go
+++ b/spatial/r3/mat_test.go
@@ -234,8 +234,9 @@ func BenchmarkQuat(b *testing.B) {
 
 var scalarFields = []struct {
 	// This is the scalar field function.
-	field   func(p Vec) float64
-	hessian func(p Vec) *Mat
+	field    func(p Vec) float64
+	gradient func(p Vec) Vec
+	hessian  func(p Vec) *Mat
 }{
 	{
 		field: func(p Vec) float64 {
@@ -243,6 +244,9 @@ var scalarFields = []struct {
 			z2 := p.Z * p.Z
 			return 4*p.X*p.X*p.X + 5*p.Y*p.Y + 3*z2*z2
 		},
+		gradient: func(p Vec) Vec {
+			return Vec{X: 12 * p.X * p.X, Y: 10 * p.Y, Z: 12 * p.Z * p.Z * p.Z}
+		},
 		hessian: func(p Vec) *Mat {
 			return NewMat([]float64{
 				24 * p.X, 0, 0,
@@ -257,6 +261,13 @@ var scalarFields = []struct {
 			y2 := p.Y * p.Y
 			return math.Cos(p.X) * math.Sin(p.Z) * y2 * y2
 		},
+		gradient: func(p Vec) Vec {
+			y3 := p.Y * p.Y * p.Y
+			y4 := p.Y * y3
+			sx, cx := math.Sincos(p.X)
+			sz, cz := math.Sincos(p.Z)
+			return Vec{X: -y4 * sx * sz, Y: 4 * y3 * cx * sz, Z: y4 * cx * cz}
+		},
 		hessian: func(p Vec) *Mat {
 			y3 := p.Y * p.Y * p.Y
 			y4 := y3 * p.Y
diff --git a/spatial/r3/vector.go b/spatial/r3/vector.go
index 29eb8342..e8e69bd2 100644
--- a/spatial/r3/vector.go
+++ b/spatial/r3/vector.go
@@ -97,6 +97,19 @@ func Divergence(p, step Vec, field func(Vec) Vec) float64 {
 	return 0.5 * (divx + divy + divz)
 }
 
+// Gradient returns the gradient of the scalar field at the point p,
+// approximated using finite differences with the given step sizes.
+func Gradient(p, step Vec, field func(Vec) float64) Vec {
+	dx := Vec{X: step.X}
+	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)),
+	}
+}
+
 // minElem return a vector with the minimum components of two vectors.
 func minElem(a, b Vec) Vec {
 	return Vec{
diff --git a/spatial/r3/vector_test.go b/spatial/r3/vector_test.go
index 8c9df4fc..7b6ce2e9 100644
--- a/spatial/r3/vector_test.go
+++ b/spatial/r3/vector_test.go
@@ -304,6 +304,25 @@ func TestDivergence(t *testing.T) {
 	}
 }
 
+func TestGradient(t *testing.T) {
+	const (
+		tol = 1e-12
+		h   = 1e-4
+	)
+	step := Vec{X: h, Y: h, Z: h}
+	rnd := rand.New(rand.NewSource(1))
+	for _, test := range scalarFields {
+		for i := 0; i < 30; i++ {
+			p := randomVec(rnd)
+			got := Gradient(p, step, test.field)
+			want := test.gradient(p)
+			if vecApproxEqual(got, want, tol) {
+				t.Errorf("result out of tolerance. got %v, want %v", got, want)
+			}
+		}
+	}
+}
+
 func vecDense(v Vec) *mat.VecDense {
 	return mat.NewVecDense(3, []float64{v.X, v.Y, v.Z})
 }