diff --git a/spatial/r3/vector.go b/spatial/r3/vector.go
index bcd0f881..9d4ab0dd 100644
--- a/spatial/r3/vector.go
+++ b/spatial/r3/vector.go
@@ -83,6 +83,18 @@ func Cos(p, q Vec) float64 {
 	return Dot(p, q) / (Norm(p) * Norm(q))
 }
 
+// Divergence returns the divergence of the vector field at the point p,
+// approximated using finite differences with the given step sizes.
+func Divergence(p, step Vec, field func(Vec) Vec) float64 {
+	sx := Vec{X: step.X}
+	divx := (field(Add(p, sx)).X - field(Sub(p, sx)).X) / step.X
+	sy := Vec{Y: step.Y}
+	divy := (field(Add(p, sy)).Y - field(Sub(p, sy)).Y) / step.Y
+	sz := Vec{Z: step.Z}
+	divz := (field(Add(p, sz)).Z - field(Sub(p, sz)).Z) / step.Z
+	return 0.5 * (divx + divy + divz)
+}
+
 // 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 ad5376b3..8c9df4fc 100644
--- a/spatial/r3/vector_test.go
+++ b/spatial/r3/vector_test.go
@@ -8,6 +8,8 @@ import (
 	"math"
 	"testing"
 
+	"golang.org/x/exp/rand"
+
 	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/mat"
 )
@@ -254,6 +256,54 @@ func TestRotate(t *testing.T) {
 	}
 }
 
+var vectorFields = []struct {
+	field      func(Vec) Vec
+	divergence func(Vec) float64
+}{
+	{
+		field: func(v Vec) Vec {
+			return Vec{X: v.X * v.Y * v.Z, Y: v.Y * v.Z, Z: v.Z * v.X}
+		},
+		divergence: func(v Vec) float64 {
+			return v.X + v.Y*v.Z + v.Z
+		},
+	},
+	{
+		field: func(v Vec) Vec {
+			sx := math.Sin(v.X)
+			sy, cy := math.Sincos(v.Y)
+			return Vec{
+				X: v.X * v.Y * v.Z * cy,
+				Y: v.Y*v.Z + sx,
+				Z: v.Z * v.X / sy,
+			}
+		},
+		divergence: func(v Vec) float64 {
+			sy, cy := math.Sincos(v.Y)
+			return v.X/sy + v.Y*v.Z*cy + v.Z
+		},
+	},
+}
+
+func TestDivergence(t *testing.T) {
+	const (
+		tol = 1e-10
+		h   = 1e-2
+	)
+	step := Vec{X: h, Y: h, Z: h}
+	rnd := rand.New(rand.NewSource(1))
+	for _, test := range vectorFields {
+		for i := 0; i < 30; i++ {
+			p := randomVec(rnd)
+			got := Divergence(p, step, test.field)
+			want := test.divergence(p)
+			if math.Abs(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})
 }