xs : [x,y,z];

beta_us : makelist(beta_u[i], i, 1, 3);

depends(beta_u, xs);

delta : ident(3);

Fs : makelist(
    apply(matrix, append([
        [-beta_u[j], 0, 0, 0, 0]
    ], makelist(
        [-diff(beta_u[j], xs[i]), -beta_u[j] * delta[i][1], -beta_u[j] * delta[i][2], -beta_u[j] * delta[i][3], -alpha * delta[i][j]]
    , i, 1, 3), 
            [
                [0, -alpha * gamma[j,1], -alpha * gamma[j,2], -alpha * gamma[j,3], -beta_u[j]]
            ])
    ),j,1,3);

results : [
    eigenvectors(Fs[1]),
    eigenvectors(Fs[2]),
    eigenvectors(Fs[3])
];

load("diag");

Lambdas : makelist(diag(apply(append, makelist(makelist(results[k][1][1][j], i, 1, results[k][1][2][j]),j,1,length(results[k][1][2])))),k,1,3);

results[3][2];