C1 : determinant ( matrix([1, 1, 1], [-1,k,-k], [-exp(-L), k*exp(k*L), -k*exp(-k*L)]) );
C2 : - determinant ( matrix([1, 1, 1], [1,k,-k], [exp(L), k*exp(k*L), -k*exp(-k*L)]) );
C3 : determinant ( matrix([1, 1, 1], [1,-1,-k], [exp(L), -exp(-L), -k*exp(-k*L)]) );
C4 : - determinant ( matrix([1, 1, 1], [1,-1,k], [exp(L), -exp(-L), k*exp(k*L)]) );

C:matrix([C1], [C2], [C3], [C4]);

D1 : determinant ( matrix([1, 1, 1], [-1,rho,-rho], [-exp(-L), rho*exp(rho*L), -rho*exp(-rho*L)]) );
D2 : - determinant ( matrix([1, 1, 1], [1,rho,-rho], [exp(L), rho*exp(rho*L), -rho*exp(-rho*L)]) );
D3 : determinant ( matrix([1, 1, 1], [1,-1,-rho], [exp(L), -exp(-L), -rho*exp(-rho*L)]) );
D4 : - determinant ( matrix([1, 1, 1], [1,-1,rho], [exp(L), -exp(-L), rho*exp(rho*L)]) );

D:matrix([D1], [D2], [D3], [D4]);

I(a,l) := a*((-1)^l *exp(a*L)-1)/(a^2 + l^2*%pi^2/L^2);

Q : matrix(
[0, -2*I(0,l),(rho-1)*I(rho+1,l),-(rho+1)*I(1-rho,l)], 
[2*I(0,l),0, (rho+1)*I(rho -1,l), (1-rho)*I(-1 - rho,l) ],
[(1-k)*I(k+1,l), - (k+1)*I(k-1,l), (rho- k)*I(k+rho,l), - (k+rho)*I(k- rho,l)],
[(k+1)*I(-k+1,l), (k-1)*I(-k-1,l) , (k+rho)*I(-k+rho,l), (k-rho)*I(-k-rho,l) ]
);

transpose(D).Q.C;