from sympy import * x, y, z, t = symbols('x y z t') def algo(X,Y) : n = len(X) P = 0 for j in range(1,int(n + 1)): L = 1 for j in range(1,int(n + 1)): if ((not (i == j))) : L = L * ((x-X[i-1])/(X[j-1]-X[i-1])) P = P + Y[i-1]*L return P print(algo([1,-3,8],[6,-8,-4]))