I'm attempting to find the eigenvectors of a matrix without using the eigenvector function.

The matrix in question is a covariance matrix:

XCov:=Matrix([[4048/5, -817/5, -122/5], [-817/5, 921/10, -1999/10], [-122/5, -1999/10, 8341/10]]);

I've already found the eigenvalues by solving for lambda:

 det := Determinant(XCov-lambda*IdentityMatrix(3));
  lambda := solve(det=0.0, lambda);

(Yes I'm reusing the eigenvalue variable for the set of eigenvalues once they've been found😏)

Anyway, I've now set up the first eigenvector I want to find as:
e1 := Vector([e11,e12,e13]);

Now, the equation to find this first eigenvector is XCov . e1 = lambda[1] . e1
I first tried putting whats on the left in a variable called eigscale(what the vector is translated to by the matrix):

eigscale := Multiply(XCov,e1);
Which returns a vector:
eigscale = [(4048/5)*e11-(817/5)*e12-(122/5)*e13,
                  -(817/5)*e11+(921/10)*e12-(1999/10)*e13,
                  -(122/5)*e11-(1999/10)*e12+(8341/10)*e13]

Each component of this vector must equate to the corresponding component in the right vector:

lambda[1]*e1 = [7.943520930*e11, 7.943520930*e12, 7.943520930*e13]

At first I tried setting these vectors equal to each other and using a solve but of course it didnt like the equations being in a vector format. So I then seperated out each equation and gave the solve function a system of equations as it expects:

solve(eigscale[1] = lambda[1]*e1[1], eigscale[2] = lambda[1]*e1[2], eigscale[3] = lambda[1]*e1[3], [e11,e12,e13]);

But again, solve fails to solve them. The reason this time(I believe) is because it can't find an exact value for e11, e12 & e13.
When solving for an eigenvector we get
e11 = e11,
e12 = Ae11,
e13 = Be11 + Ce12

I was wondering if there was a way to do a partial solve to find the components in terms of each other?

Failing that, I'm aware I can do it manually through row operations but I believe that would require changing the format so that each equation is a component of a single vector:
eigsolve := Vector([eigscale[1] = lambda[1]*e1[1], eigscale[2] = lambda[1]*e1[2], eigscale[3] = lambda[1]*e1[3]]);

Since row operations cannot be performed on a equation of vectors (again, I believe).

Help appreciated!

So I was trying to create a shorthand for creating a plot of multiple arrows, with the arrow colour dependent on the magnitude of the vector.
I currently have a set of vectors, v, I want to display, and v[4] is the largest.