#computing the truncated tau function for the (3,6)-soliton in Example 8
restart:
k := 3:
n := 6:
A := Matrix([[1, 1, 0, 0, 0, 0], [0, 0, 1, 1, 0, 0], [0, 0, 0, 0, 1, 1]]):


K := Array(1 ..n):
for i to n do
    K[i] := kappa(i)
end do:


N := 20:
sizeBasis := 15:
f := Array(1 ..sizeBasis):

#compute a finite number of the functions which give a basis for the point in the Sato Grassmannian 
i := 1:
for j to k do
    for i to n do 
        f[j] := f[j] + A[j][i]/(-z*K[i] + 1)
end do:
    f[j] := expand(factor(f[j]*z^(-k + 1)*z^(sizeBasis - 1))):
end do:


for j to sizeBasis - k do
    f[j + k] := z^(-j - k + 1)*z^(sizeBasis - 1)
end do:

#build and display the matrix
Cmatrix := Array(1 ..sizeBasis):

for i to sizeBasis do
    Cmatrix[i] := mtaylor(f[i], z, N + i)
end do:

xi := (i, j) -> factor(coeff(Cmatrix[sizeBasis - j + 1], z, i - 1)):
A := Matrix(N, sizeBasis, xi);


#compute the elementary Schur polynomials in three variables up to partitions of partLength 

partLength := 8:
p = Array(1 ..partLength):
for i to partLength do
    p[i] := 0
end do:
n1 := 0: n2 := 0: n3 := 0:
for i to partLength do
    for n1 from 0 to i do
        for n2 from 0 to i do
            for n3 from 0 to i do
                if n1 + 2*n2 + 3*n3 = i then
                   p[i] := p[i] + X^n1*Y^n2*T^n3/(n1!*n2!*n3!)
                end if: 
            end do:
        end do:
    end do:
end do:


with(combinat, partition):with(ListTools):
partitions := FlattenOnce([seq(partition(j), j = 1 ..partLength)]):
h := numelems(partitions):

#compute the Maya diagrams corresponding to the partitions
Maya = Array(1 .. h):
for i to h do
    m := numelems(partitions[I]):
    Maya[i] := Array(1 .. m):
    for j to m do
        Maya[i][j] := partitions[i][m - j + 1] - j + sizeBasis + 1:
    end do:
end do:
tau := 0:

#compute the coordinates, the Schur polynomial for the non zero coordinates and the tau function
Coordinates = Array(1 .. h):
with(LinearAlgebra): with(ArrayTools): with(linalg):
for i to h do
    m := numelems(Maya[i]):
    increasingMaya := convert(FlipDimension(Maya[i], 2), list):
    if Maya[i][1] < N then
        Coordinates[i] := Determinant(SubMatrix(A, increasingMaya, [sizeBasis - m + 1 ..sizeBasis])):
    else
        Coordinates[i] := 0:
    end if:
    if Coordinates[i] <> 0 then
        l := Array(1 ..m):
        v := Vector(1 ..m):
        for k to m do
            l[k] := partitions[i][k] + k - 1:
            v[k] := p[l[k]]:
        end do:
        Wr := wronskian(v, X):
        schur := det(Wr):
        tau := schur*Coordinates[i] + tau:
    end if:
end do:
tau;