diff --git a/scripts/spline.py b/scripts/spline.py
new file mode 100644
index 0000000..41c5465
--- /dev/null
+++ b/scripts/spline.py
@@ -0,0 +1,93 @@
+import numpy as np
+from math import sqrt
+
+def cubic_interp1d(x0, x, y):
+    """
+    Interpolate a 1-D function using cubic splines.
+      x0 : a float or an 1d-array
+      x : (N,) array_like
+          A 1-D array of real/complex values.
+      y : (N,) array_like
+          A 1-D array of real values. The length of y along the
+          interpolation axis must be equal to the length of x.
+
+    Implement a trick to generate at first step the cholesky matrice L of
+    the tridiagonal matrice A (thus L is a bidiagonal matrice that
+    can be solved in two distinct loops).
+
+    additional ref: www.math.uh.edu/~jingqiu/math4364/spline.pdf 
+    """
+    x = np.asfarray(x)
+    y = np.asfarray(y)
+
+    # remove non finite values
+    # indexes = np.isfinite(x)
+    # x = x[indexes]
+    # y = y[indexes]
+
+    # check if sorted
+    if np.any(np.diff(x) < 0):
+        indexes = np.argsort(x)
+        x = x[indexes]
+        y = y[indexes]
+
+    size = len(x)
+
+    xdiff = np.diff(x)
+    ydiff = np.diff(y)
+
+    # allocate buffer matrices
+    Li = np.empty(size)
+    Li_1 = np.empty(size-1)
+    z = np.empty(size)
+
+    # fill diagonals Li and Li-1 and solve [L][y] = [B]
+    Li[0] = sqrt(2*xdiff[0])
+    Li_1[0] = 0.0
+    B0 = 0.0 # natural boundary
+    z[0] = B0 / Li[0]
+
+    for i in range(1, size-1, 1):
+        Li_1[i] = xdiff[i-1] / Li[i-1]
+        Li[i] = sqrt(2*(xdiff[i-1]+xdiff[i]) - Li_1[i-1] * Li_1[i-1])
+        Bi = 6*(ydiff[i]/xdiff[i] - ydiff[i-1]/xdiff[i-1])
+        z[i] = (Bi - Li_1[i-1]*z[i-1])/Li[i]
+
+    i = size - 1
+    Li_1[i-1] = xdiff[-1] / Li[i-1]
+    Li[i] = sqrt(2*xdiff[-1] - Li_1[i-1] * Li_1[i-1])
+    Bi = 0.0 # natural boundary
+    z[i] = (Bi - Li_1[i-1]*z[i-1])/Li[i]
+
+    # solve [L.T][x] = [y]
+    i = size-1
+    z[i] = z[i] / Li[i]
+    for i in range(size-2, -1, -1):
+        z[i] = (z[i] - Li_1[i-1]*z[i+1])/Li[i]
+
+    # find index
+    index = x.searchsorted(x0)
+    np.clip(index, 1, size-1, index)
+
+    xi1, xi0 = x[index], x[index-1]
+    yi1, yi0 = y[index], y[index-1]
+    zi1, zi0 = z[index], z[index-1]
+    hi1 = xi1 - xi0
+
+    # calculate cubic
+    f0 = zi0/(6*hi1)*(xi1-x0)**3 + \
+         zi1/(6*hi1)*(x0-xi0)**3 + \
+         (yi1/hi1 - zi1*hi1/6)*(x0-xi0) + \
+         (yi0/hi1 - zi0*hi1/6)*(xi1-x0)
+    return f0
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    x = np.linspace(0, 10, 11)
+    y = np.sin(x)
+    plt.scatter(x, y)
+
+    x_new = np.linspace(0, 10, 201)
+    plt.plot(x_new, cubic_interp1d(x_new, x, y))
+
+    plt.show()