from numpy import *
import matplotlib.pyplot as plt

x = array([[0],[1],[2]])
f = array([[1],[3],[2]])

n = shape(x)[0]-1
V = x**(arange(n+1))
c = linalg.solve(V,f)

x_plot = linspace(x[0,0]-0.1, x[-1,0]+0.1, 100).reshape((100,1))
x_plot_powers = x_plot**(arange(n+1));
plt.plot(x_plot, x_plot_powers@c)
plt.show()

n = shape(x)[0] - 1
A = zeros((n+1,n+1))
A[:,0] = 1;
for c in range(1,n+1):
    A[c:,c] = A[c:,c-1]*(x[c:,0]-x[c-1,0])
    
c = zeros((n+1,1));
for k in range(n+1):
    c[k,0] = (f[k,0] - A[k,:k]@c[:k,0])/A[k,k]

x_plot = linspace(x[0,0]-0.1,x[-1,0]+0.1,100).reshape((100,1))
x_plot_powers = ones((shape(x_plot)[0], n+1))
for k in range(1,n+1):
    x_plot_powers[:,k] = x_plot_powers[:,k-1]*(x_plot[:,0]-x[k-1,0])

plt.plot(x_plot, x_plot_powers@c)
plt.show()