diff --git a/govector_test.go b/govector_test.go
index ea08469..077b64c 100644
--- a/govector_test.go
+++ b/govector_test.go
@@ -562,3 +562,89 @@ func TestPointInPolygon(t *testing.T) {
 		t.Error("point should have been in polygon")
 	}
 }
+
+func TestExtrapolate(t *testing.T) {
+	_, err := Extrapolate([]Point2d{})
+	if err == nil {
+		t.Error("should require len(P) == 3")
+	}
+
+	tests := []struct {
+		points   []Point2d
+		expected Point2d
+	}{
+		{
+			points: []Point2d{
+				{0, 0},
+				{1, 1},
+				{2, 2},
+			},
+			expected: Point2d{3, 3},
+		},
+		{
+			points: []Point2d{
+				{0, 0},
+				{1, 1},
+				{2, 1},
+			},
+			expected: Point2d{3, 0},
+		},
+		{
+			points: []Point2d{
+				{0, 0},
+				{1, 0},
+				{2, 0},
+			},
+			expected: Point2d{3, 0},
+		},
+		{
+			points: []Point2d{
+				{0, 0},
+				{0, 1},
+				{0, 2},
+			},
+			expected: Point2d{0, 3},
+		},
+		{
+			points: []Point2d{
+				{0, 0},
+				{0, 0},
+				{0, 0},
+			},
+			expected: Point2d{0, 0},
+		},
+		{
+			points: []Point2d{
+				{42, 42},
+				{42, 42},
+				{42, 42},
+			},
+			expected: Point2d{42, 42},
+		},
+		{
+			points: []Point2d{
+				{104, 288},
+				{110, 270},
+				{120, 250},
+			},
+			expected: Point2d{134, 228},
+		},
+		{
+			points: []Point2d{
+				{0, 0},
+				{5, 3},
+				{7, 3},
+			},
+			expected: Point2d{6, 0},
+		},
+	}
+	for _, test := range tests {
+		p, err := Extrapolate(test.points)
+		if err != nil {
+			t.Errorf("yeah, this should've been nil: %s", err)
+		}
+		if *p != test.expected {
+			t.Errorf("wrong: expected: %+v, actual: %+v", test.expected, p)
+		}
+	}
+}
diff --git a/interp.go b/interp.go
new file mode 100644
index 0000000..0322983
--- /dev/null
+++ b/interp.go
@@ -0,0 +1,23 @@
+package govector
+
+import (
+	"errors"
+)
+
+// explicitly implements Neville's algorithm for three points at very
+// particular parameter values.
+// TODO: update this to take an arbitrary parameter value, and arbitrary slice
+// of Point2d
+func Extrapolate(P []Point2d) (*Point2d, error) {
+	if len(P) != 3 {
+		return nil, errors.New("only works for len(P) == 3")
+	}
+
+	// P_{01} = -2 * P_0 + 3 * P_1
+	p01 := P[0].ToVector().Scale(-2.0).ToPoint().Add(P[1].ToVector().Scale(3.0))
+	// P_{12} = -1 * P_1 + 2 * P_2
+	p12 := P[1].ToVector().Scale(-1.0).ToPoint().Add(P[2].ToVector().Scale(2.0))
+	// P_{012} = 1/2 * P_{01} + 3/2 * P_{12}
+	p012 := p01.ToVector().Scale(-1.0 / 2.0).ToPoint().Add(p12.ToVector().Scale(3.0 / 2.0))
+	return &p012, nil
+}