load("eigen");

load("diag");

eigenvalueList(results) := apply(append,makelist(makelist(results[1][1][j],i,1,results[1][2][j]),j,1,length(results[1][1])));

rightEigenvectorMatrix(results) := transpose(apply(matrix, apply(append, results[2])));

assume(%alpha>0, %gamma[xx]>0, f>0);

F : matrix(
    [0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0],
    [%alpha^2 * K[xx] / %gamma[xx] * (f + %alpha * df),
        -%alpha * K[xx] * f / %gamma[xx]^2,
        0, 0,
        %alpha * f / %gamma[xx]],
    [%alpha^2 * K[xx], 0, 0, 0, %alpha],
    [%alpha^2 * a[x], 0, %alpha, 0, 0]
);

results : eigenvectors(F);

eigenvalues : eigenvalueList(results);

eigenvalueMatrix : diag(eigenvalues);

rightEigenvectorMatrix(results)$
rightEigenvectors : %;

transpose(%)$
makelist(transpose(matrix(%[i])),i,1,4)$
apply(orthogonal_complement, %)$
addcol(rightEigenvectors, args(%)[1])$
/*addcol(rightEigenvectors, transpose(matrix([0,1,0,0,0])))$*/
rightEigenvectors : %;

leftEigenvectors : ratsimp(invert(rightEigenvectors));

rightEigenvectors . eigenvalueMatrix . leftEigenvectors - F$
ratsimp(%);

leftEigenvectorsProvided : matrix(
    [0, 0, 1/sqrt(%gamma[xx]), 0, sqrt(f)/%gamma[xx]],
    [0, 0, -1/sqrt(%gamma[xx]), 0, sqrt(f)/%gamma[xx]],
    [1, 0, 0, 0, 0],
    [0, 1, 0, 0, 0],
    [0, 0, 1, -f/%gamma[xx], 0]
);

rightEigenvectorsProvided : ratsimp(invert(leftEigenvectorsProvided));

eigenvalueMatrixProvided : diag([%alpha * sqrt(f/%gamma[xx]), -%alpha * sqrt(f/%gamma[xx]), 0, 0, 0]);

rightEigenvectorsProvided . eigenvalueMatrixProvided . leftEigenvectorsProvided;