diff --git a/python/NevilleInterpolation.py b/python/NevilleInterpolation.py
new file mode 100644
index 0000000..c71f3fb
--- /dev/null
+++ b/python/NevilleInterpolation.py
@@ -0,0 +1,29 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from NevillePoint import nevillepoint
+
+
+# [USES] interpolations/NevillePoint
+# Function that calculates the Neville interpolation polynomial for a vector
+# of points (a, b) with a:xs, b:ys in "elements" points between the first of xs
+# and the last of xs; it also plots the interpolation
+def nevilleinterpolation(xs, ys, elements=100):
+    n = len(xs)
+    xss = np.linspace(xs[0], xs[n - 1], elements)
+    # create an empty vector to store the values of the interpolation in the
+    # xss' points
+    yss = []
+
+    for i in range(0, len(xss)):
+        # Calculate the approximation of the function in xss(i)
+        yss.append(nevillepoint(0, n - 1, xss[i], xs, ys))
+
+
+
+    print(yss)
+    # plot the interpolated function
+    plt.plot(xss, yss)
+    # plot the initial points as circles
+    plt.plot(xs, ys, 'o', "linewidth", 3)
+    plt.show()
+
diff --git a/python/NevillePoint.py b/python/NevillePoint.py
new file mode 100644
index 0000000..40f1838
--- /dev/null
+++ b/python/NevillePoint.py
@@ -0,0 +1,22 @@
+# Function that calculates the value of the Neville interpolation polynomial
+# P_ij(x) for a vector of points (a, b) with a:xs, b:xs in the abscissa x;
+# [NOTE] It doesn't use explicit caching (just the internal Octave caching
+# system)
+
+
+def nevillepoint(i, j, x, xs, ys):
+    # Base case: P_ii(x) = f(x_i)
+    # ys is 1 - indexed, so we have to add one to the index
+    if i == j:
+       y = ys[i]
+       return y
+
+    # x_j - x_i
+    delta = xs[j] - xs[i]
+    # P(i, j - 1, x)
+    pij_ = (xs[j] - x) * nevillepoint(i, j - 1, x, xs, ys)
+    # P(i + 1, j, x)
+    pi_j = (x - xs[i]) * nevillepoint(i + 1, j, x, xs, ys)
+    # P(i, j, x)
+    y = (pij_ + pi_j) / delta
+    return y