import numpy as np
from scipy.integrate import quad as integrator

def approximate(f, n, x0, x1):

    A = np.zeros((n,n))
    b = np.zeros((n,1))

    for i in range(n):
        b[i] = integrator(lambda x: f(x)*x**i, x0, x1)[0] #slide 74, Ch 5.
        for j in range(n):
            A[i,j] = integrator(lambda x: (x**i)*(x**j), x0, x1)[0]


    p = np.linalg.solve(A, b)

    return p

def f(x):
    return np.sin(x)
   

def eval(p, x):

    y = 0*x; 
    for i in range(len(p)):
	y = y + p[i]*(x**i)

    return y


x0 = 0
x1 = 2*np.pi
p = approximate(f, 4, x0, x1)

print p



x = np.linspace(x0, x1, 1001)
y = eval(p, x)


import pylab
pylab.clf()
pylab.plot(x, y, x,f(x))
pylab.show()