In [1]: 
from __future__ import division

import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as pt

Mini-Introduction to sympy

In [2]: 
import sympy as sym

# Enable "pretty-printing" in IPython
sym.init_printing()
In [9]: 
x = sym.Symbol("x")

myexpr = (x**2-3)**2
myexpr
Out[9]:
$$\left(x^{2} - 3\right)^{2}$$
In [11]: 
myexpr.expand()
Out[11]:
$$x^{4} - 6 x^{2} + 9$$
In [10]: 
sym.integrate(myexpr, x)
Out[10]:
$$\frac{x^{5}}{5} - 2 x^{3} + 9 x$$
In [16]: 
sym.integrate(myexpr, (x, -1, 1))
Out[16]:
$$\frac{72}{5}$$

Orthogonal polynomials

In [22]: 
def inner_product(f1, f2):
    return sym.integrate(f1*f2, (x, -1, 1))
In [23]: 
inner_product(1, 1)
Out[23]:
$$2$$
In [24]: 
inner_product(1, x)
Out[24]:
$$0$$
In [63]: 
#basis = [1, x, x**2, x**3]
basis = [1, x, x**2, x**3, x**4, x**5]
In [64]: 
# Now for plain old Gram-Schmidt
orth_basis = []

for q in basis:
    for prev_q in orth_basis:
        q = q - inner_product(prev_q, q)*prev_q
    q = q / sym.sqrt(inner_product(q, q))
    
    orth_basis.append(q)
In [65]: 
orth_basis
Out[65]:
$$\begin{bmatrix}\frac{\sqrt{2}}{2}, & \frac{\sqrt{6} x}{2}, & \frac{3 \sqrt{10}}{4} \left(x^{2} - \frac{1}{3}\right), & \frac{5 \sqrt{14}}{4} \left(x^{3} - \frac{3 x}{5}\right), & \frac{105 \sqrt{2}}{16} \left(x^{4} - \frac{6 x^{2}}{7} + \frac{3}{35}\right), & \frac{63 \sqrt{22}}{16} \left(x^{5} - \frac{10 x^{3}}{9} + \frac{5 x}{21}\right)\end{bmatrix}$$

(Note: normalization different than in the book.)

In [66]: 
mesh = np.linspace(-1, 1, 100)

pt.figure(figsize=(8,8))
for f in orth_basis:
    f = sym.lambdify(x, f)
    pt.plot(mesh, [f(xi) for xi in mesh])
In []: