v : transpose(matrix([v0, v1, v2, v3]));

A_x : matrix([0, -c * gxx, -c * gxy, -c * gxz], [-c,0,0,0], [0,0,0,0], [0,0,0,0]);

A_y : matrix([0, -c * gxy, -c * gyy, -c * gyz], [0,0,0,0], [-c,0,0,0], [0,0,0,0]);

A_z : matrix([0, -c * gxz, -c * gyz, -c * gzz], [0,0,0,0], [0,0,0,0], [-c,0,0,0]);

results : eigenvectors(A_x);

load(diag);

P : matrix(
    [1, 0, 0, 0],
    [0, 0, 1, 0],
    [0, 0, 0, 1],
    [0, 1, 0, 0]
);
IP : invert(P);

Lambda : P . diag([results[1][1][1], results[1][1][2], results[1][1][3], results[1][1][3]]) . IP;

R : transpose(matrix(results[2][1][1], results[2][2][1], results[2][3][1], results[2][3][2])) . diag([sqrt(gxx), sqrt(gxx), gxz, gxz]) . IP;

L : ratsimp(invert(R));

R . Lambda . L$
ratsimp(%);

/* From Wolfram Alpha */;

P : matrix(
    [0, 0, 1, 0],
    [1, 0, 0, 0],
    [0, 1, 0, 0],
    [0, 0, 0, 1]
);
IP : invert(P);

Lambda_x : P . diag([0, 0, -c * sqrt(gxx), c * sqrt(gxx)]) . IP;

R_x : transpose(matrix(
    [0, -gxy, gxx, 0],
    [0, -gxz, 0, gxx],
    [sqrt(gxx), 1, 0, 0],
    [-sqrt(gxx), 1, 0, 0]
    )) . IP;

L_x : ratsimp(invert(R_x));

R_x . Lambda_x . L_x$ ratsimp(%);

Lambda_y : P . diag([0, 0, -c * sqrt(gyy), c * sqrt(gyy)]) . IP;

R_y : transpose(matrix(
    [0, 0, -gyz, gyy],
    [0, gyy, -gxy, 0],
    [sqrt(gyy), 0, 1, 0],
    [-sqrt(gyy), 0, 1, 0]
)) . IP;

L_y : ratsimp(invert(R_y));

R_y . Lambda_y . L_y$ ratsimp(%);

Lambda_z : diag([-c * sqrt(gzz), 0, 0, c * sqrt(gzz)]);

R_z : transpose(matrix(
    [0, gzz, 0, -gxz],
    [0, 0, gzz, -gyz],
    [sqrt(gzz), 0, 0, 1],
    [-sqrt(gzz), 0, 0, 1]
)) . IP;

L_z : ratsimp(invert(R_z));

R_z . Lambda_z . L_z$ ratsimp(%);

L_x . v$ ratsimp(%);

L_y . v$ ratsimp(%);

L_z . v$ ratsimp(%);

R_x . v$ ratsimp(%);

R_y . v$ ratsimp(%);

R_z . v$ ratsimp(%);

A_x . R_x = R_x . Lambda_x$ ratsimp(%); is(%);

A_y . R_y = R_y . Lambda_y$ ratsimp(%); is(%);

A_z . R_z = R_z . Lambda_z$ ratsimp(%); is(%);