SymPy

I've been using Wolfram Alpha for my integration, derivation and algebraic manipulation needs, but I tried SymPy and I have to say that it's pretty neat. I used it to implement the Gram-Schmidt Process to obtain a set of orthonormal basis functions from a set of non-orthogonal, non-normalized and linearly independent polynomial functions and SymPy has done the job beautifully. If you already know a little Python, the learning curve is quite gentle.

from __future__ import division
from sympy import *
from prettytable import PrettyTable
import functools
import string

x=symbols('x')
init_printing(use_unicode=True)

def print_basis(basis):
    """
    pretty print basis functions
    """
    table = PrettyTable(["e%d" %i for i in range(len(basis))])
    table.padding_width = 1
    table.add_row(map(pretty, basis))
    print table

def print_projection_coefficients(coefficients, significant_digits):
    """
    pretty print projection coefficients
    """
    table = PrettyTable(string.ascii_uppercase[:len(coefficients)])
    table.padding_width = 1
    table.add_row(map(pretty, coefficients))
    table.add_row(map(functools.partial(N, n=significant_digits), coefficients))
    print table

def normalize(f, inner_product):
    """
    normalize input function using input inner product operation
    """
    return f/sqrt(inner_product(f, f))

def project(f, basis):
    """
    project input function in space generated by input basis
    """
    projection = 0
    coefficients = []
    for b in basis:
        coefficient = inner_product(f, b)
        coefficients.append(coefficient)
        projection = projection + coefficient * b
    return projection, coefficients

def gram_schmidt(v_basis, inner_product):
    """
    orthonormalize input basis vectors using gram-schmidt algorithm
    """
    e = normalize(v_basis[0], inner_product)
    e_basis = [e]
    for v in v_basis[1:]:
        f = v
        for e in e_basis:
            f = f - inner_product(e, v) * e
        e_basis.append(normalize(f, inner_product))
    return e_basis

# define non orthonormal basis functions here
v_basis = [1, x, x**2, x**3, x**4]

# define inner_product here
inner_product=lambda f1, f2: integrate(f1*f2,(x, -pi, pi))

# orthonormalize basis using gram-schmidt algorithm
e_basis = gram_schmidt(v_basis, inner_product)
print("orthonormalized basis functions:")
print_basis(e_basis)

# get projection of cos(x) in orthonormal basis
significant_digits=6
approximation, coefficients = project(cos(x), e_basis)
print("projection coefficients of cos(x):")
print_projection_coefficients(coefficients, significant_digits)
print("projection of cos(x):")
pprint(approximation.evalf(significant_digits))

SymPy formats the output nicely in unicode. I used PrettyTable to print the output nicely in table form. I'm very pleased with the result. Together with NumPy and Matplotlib, it makes the perfect combo for scientific programming.

orthonormalized basis functions:
+---------+---------+--------------------+------------------------+--------------------------------------+
|    e0   |    e1   |         e2         |           e3           |                  e4                  |
+---------+---------+--------------------+------------------------+--------------------------------------+
|    ___  |   ___   |          ⎛      2⎞ |          ⎛        2  ⎞ |           ⎛          ⎛      2⎞     ⎞ |
|  ╲╱ 2   | ╲╱ 6 ⋅x |     ____ ⎜ 2   π ⎟ |     ____ ⎜ 3   3⋅π ⋅x⎟ |           ⎜        2 ⎜ 2   π ⎟     ⎟ |
| ─────── | ─────── | 3⋅╲╱ 10 ⋅⎜x  - ──⎟ | 5⋅╲╱ 14 ⋅⎜x  - ──────⎟ |           ⎜     6⋅π ⋅⎜x  - ──⎟    4⎟ |
|     ___ |     3/2 |          ⎝     3 ⎠ |          ⎝       5   ⎠ |       ___ ⎜ 4        ⎝     3 ⎠   π ⎟ |
| 2⋅╲╱ π  |  2⋅π    | ────────────────── | ────────────────────── | 105⋅╲╱ 2 ⋅⎜x  - ────────────── - ──⎟ |
|         |         |          5/2       |            7/2         |           ⎝           7          5 ⎠ |
|         |         |       4⋅π          |         4⋅π            | ──────────────────────────────────── |
|         |         |                    |                        |                   9/2                |
|         |         |                    |                        |               16⋅π                   |
+---------+---------+--------------------+------------------------+--------------------------------------+
projection coefficients of cos(x):
+---+---+------------+---+----------------------------------------------------------+
| A | B |     C      | D |                            E                             |
+---+---+------------+---+----------------------------------------------------------+
| 0 | 0 |      ____  | 0 |             ⎛            3⎞             ⎛      3       ⎞ |
|   |   | -3⋅╲╱ 10   |   |         ___ ⎜        16⋅π ⎟         ___ ⎜  16⋅π        ⎟ |
|   |   | ────────── |   |   105⋅╲╱ 2 ⋅⎜-24⋅π + ─────⎟   105⋅╲╱ 2 ⋅⎜- ───── + 24⋅π⎟ |
|   |   |     3/2    |   |             ⎝          7  ⎠             ⎝    7         ⎠ |
|   |   |    π       |   | - ───────────────────────── + ────────────────────────── |
|   |   |            |   |                9/2                         9/2           |
|   |   |            |   |            16⋅π                        16⋅π              |
| 0 | 0 |  -1.70371  | 0 |                         0.486659                         |
+---+---+------------+---+----------------------------------------------------------+
projection of cos(x):
           4             2           
0.0261599⋅x  - 0.452289⋅x  + 0.978327