In sage, there is a method `gram_schnidt`

associated to a matrix.

After the matrix is declared, the method can be immediately called.
The following defines a matrix `A`

, computes the Gram-Schnidt matrices `G.M`

as in the examples from

```
A.gram_schmidt?
```

and does a check for the relation $A=MG$, it is a possible starting point:

```
var( 'x' )
def ip(f, g):
return integrate( f * g, x, -1, 1 )
f_list = [1, x, x^2, x^3, x^4]
A = matrix( QQ, len(f_list),
[ ip( f, g ) for f in f_list for g in f_list ] )
print "A is\n%s\n" % A
G, M = A.gram_schmidt()
print "G is\n%s\n" % G
print "M is\n%s\n" % M
print "Is MG = A? %s" % bool( M*G == A )
```

Results:

```
A is
[ 2 0 2/3 0 2/5]
[ 0 2/3 0 2/5 0]
[2/3 0 2/5 0 2/7]
[ 0 2/5 0 2/7 0]
[2/5 0 2/7 0 2/9]
G is
[ 2 0 2/3 0 2/5]
[ 0 2/3 0 2/5 0]
[ -424/5439 0 1376/9065 0 248/1813]
[ 0 -12/595 0 4/119 0]
[ 64/112035 0 -128/22407 0 64/9603]
M is
[ 1 0 0 0 0]
[ 0 1 0 0 0]
[675/1813 0 1 0 0]
[ 0 75/119 0 1 0]
[425/1813 0 950/1067 0 1]
Is MG = A? True
```

**Later edit:**

I try now hard to reproduce in sagemath the following result from maxima,
obtained by copy+pasting into the maxima interpreter the code,
thanks go to danielvolinski for this sample,

```
load(eigen);
ip(f,g) := integrate(f*g, u, 0, 1);
y : gramschmidt([1, u, u^2, u^3, u^4], ip);
```

which delivers

```
2 2
2 u - 1 6 u - 6 u + 1 (2 u - 1) (10 u - 10 u + 1)
(%o14) [1, -------, --------------, ----------------------------,
2 6 20
4 3 2
70 u - 140 u + 90 u - 20 u + 1
---------------------------------]
70
```

The situation is different, because there is no *a priori* knowledge of an orthonormal basis of the vector space generated by $1,u,u^2,u^3,u^4$ with the given scalar product $\langle f,g\rangle =\int_0^1f\;g$ .

To get a similar result i prefer to program from the scratch, there is no implemented facility...:

```
var( 'u' )
def ip(f, g):
return integrate( f * g, u, 0, 1 )
n = 5
B = vector( [u^k for k in range(n)] )
T = identity_matrix( QQ, n ) # in T we will have finally the transfromation matrix
C = vector( [u^k for k in range(n)] ) # in C we will have finally the Gram-Schmidt basis
for k in range(n-1):
S = identity_matrix( QQ, n )
for j in range(k+1, n):
S[j, k] = -ip( C[j], C[k] ) / ip( C[k], C[k] )
T = S*T
C = S*C
print T
print C
```

This gives the transformation matrix, and the new basis

```
[ 1 0 0 0 0]
[ -1/2 1 0 0 0]
[ 1/6 -1 1 0 0]
[-1/20 3/5 -3/2 1 0]
[ 1/70 -2/7 9/7 -2 1]
(1, u - 1/2, u^2 - u + 1/6, u^3 - 3/2*u^2 + 3/5*u - 1/20, u^4 - 2*u^3 + 9/7*u^2 - 2/7*u + 1/70)
```

in concordance with the maxima result.