diff --git a/govector_test.go b/govector_test.go
index 6fde00b..dedf929 100644
--- a/govector_test.go
+++ b/govector_test.go
@@ -687,6 +687,22 @@ func TestNeville(t *testing.T) {
 		t        float64
 		expected Point2d
 	}{
+		{
+			points: []Point2d{
+				{0, 0},
+				{3, 3},
+				{5, 3},
+				{7, 0},
+			},
+			ts: []float64{
+				0,
+				1,
+				3,
+				4,
+			},
+			t:        2.0,
+			expected: Point2d{25.0 / 6.0, 4.0},
+		},
 		{
 			points: []Point2d{
 				{0, 0},
@@ -694,13 +710,97 @@ func TestNeville(t *testing.T) {
 				{2, 2},
 			},
 			ts: []float64{
-				0.0,
+				0,
 				1.0 / 3.0,
 				2.0 / 3.0,
 			},
 			t:        1.0,
 			expected: Point2d{3, 3},
 		},
+		{
+			points: []Point2d{
+				{0, 0},
+				{1, 1},
+				{2, 1},
+			},
+			ts: []float64{
+				0,
+				1.0 / 3.0,
+				2.0 / 3.0,
+			},
+			t:        1.0,
+			expected: Point2d{3, 0},
+		},
+		{
+			ts: []float64{
+				0,
+				1.0 / 3.0,
+				2.0 / 3.0,
+			},
+			t: 1.0,
+			points: []Point2d{
+				{0, 0},
+				{1, 0},
+				{2, 0},
+			},
+			expected: Point2d{3, 0},
+		},
+		{
+			points: []Point2d{
+				{0, 0},
+				{0, 1},
+				{0, 2},
+			},
+			ts: []float64{
+				0,
+				1.0 / 3.0,
+				2.0 / 3.0,
+			},
+			t:        1.0,
+			expected: Point2d{0, 3},
+		},
+		{
+			ts: []float64{
+				0,
+				1.0 / 3.0,
+				2.0 / 3.0,
+			},
+			t: 1.0,
+			points: []Point2d{
+				{0, 0},
+				{0, 0},
+				{0, 0},
+			},
+			expected: Point2d{0, 0},
+		},
+		{
+			points: []Point2d{
+				{42, 42},
+				{42, 42},
+				{42, 42},
+			},
+			expected: Point2d{42, 42},
+			ts: []float64{
+				0,
+				1.0 / 3.0,
+				2.0 / 3.0,
+			},
+			t: 1.0,
+		},
+		{
+			points: []Point2d{
+				{104, 288},
+				{110, 270},
+				{120, 250},
+			},
+			expected: Point2d{134, 228},
+			ts: []float64{
+				0,
+				1.0 / 3.0,
+				2.0 / 3.0,
+			},
+			t: 1.0,
+		},
 	}
 	for _, test := range tests {
 		p, err := Neville(test.points, test.ts, test.t)
@@ -708,8 +808,13 @@ func TestNeville(t *testing.T) {
 			t.Errorf("yeah, this should've been nil: %s", err)
 		}
 		if p == nil {
-		} else if *p != test.expected {
-			t.Errorf("wrong: expected: %+v, actual: %+v", test.expected, p)
+			t.Errorf("nil point?")
+		} else {
+			deltaX := math.Abs(float64(p.X - test.expected.X))
+			deltaY := math.Abs(float64(p.Y - test.expected.Y))
+			if deltaX > 1e-4 || deltaY > 1e-4 { // HOLY SHIT!!!
+				t.Errorf("wrong: expected: %+v, actual: %+v", test.expected, p)
+			}
 		}
 	}
 }
diff --git a/interp.go b/interp.go
index 40dad2e..444ad01 100644
--- a/interp.go
+++ b/interp.go
@@ -6,9 +6,7 @@ import (
 )
 
 // 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
+// particular parameter values: ts[]{0, 1/3, 2/3}, extrapolate to 1.
 func Extrapolate(P []Point2d) (*Point2d, error) {
 	if len(P) != 3 {
 		return nil, errors.New("only works for len(P) == 3")
@@ -23,6 +21,8 @@ func Extrapolate(P []Point2d) (*Point2d, error) {
 	return &p012, nil
 }
 
+// Implements Neville's Algorithm for a slice of points (P) at parameter
+// values(ts), to parameter value (t). Returns the point and an error
 func Neville(P []Point2d, ts []float64, t float64) (*Point2d, error) {
 	if len(P) != len(ts) {
 		return nil, errors.New(
@@ -33,5 +33,16 @@ func Neville(P []Point2d, ts []float64, t float64) (*Point2d, error) {
 			),
 		)
 	}
-	return nil, errors.New("smb")
+	Q := []Point2d(P)
+	order := len(P) - 1
+	for i := 0; i < order; i++ {
+		for j := 0; j < order-i; j++ {
+			var tLeft float32 = float32(ts[i+j+1] - t)
+			var tRight float32 = float32(t - ts[j])
+			Q[j] = Q[j].ToVector().Scale(tLeft).ToPoint().Add(P[j+1].ToVector().Scale(tRight))
+			Q[j] = Q[j].ToVector().Scale(float32(1.0 / (ts[i+j+1] - ts[j]))).ToPoint()
+
+		}
+	}
+	return &Q[0], nil
 }