def G_ls (xi, basis):

    M = numpy.size (xi)
    N = numpy.size (basis)
    G = numpy.zeros (N*N).reshape(N, N)

    for k in xrange (N):
        for l in xrange (N):
            G[k, l] = numpy.sum ((basis[l](xi[i]) * basis[k](xi[i])) for i in xrange (M))
    return G