Thank you for the response, Markiyan.

First of all, there should be a zero in the last row, second column of the original matrix I posted. Sorry about that mistake. However, your suggestion solved my problem and allowed me to calculate the eigenvectors. So thank you very much for that tip.

The trigonometric functions still cause some issues in my answer. Here is one of the eigenvectors I get:

u := Vector(3, {(1) = 489794*sin(theta)*cos(theta)/(497206*cos(theta)^2-245103-sqrt(7315644000*cos(theta)^4-3835202000*cos(theta)^2+60075480609)), (2) = 0, (3) = 1})

If I write eval(u[1], theta=0), I get zero as the answer. However, by a quick check I see that both the numerator and denominator are both zero for theta=0. Letting theta go to zero with the help of limit(u[1], theta=0) makes the result 'undefined' (I guess this means answer goes to infinity), as it should be. First question is: why doesn't eval() recognize that there is division by zero?

I have checked by hand that the eigenvectors are correct. Second question: is there some way to avoid this division by zero? I have tried for example to scale the eigenvectors, but that led to a zero-vector. I also tried to "invert" the eigenvector, such that it became

u := Vector(3, {(1) =1, (2) = 0, (3) = (1/489794)*(497206*cos(theta)^2-245103-sqrt(7315644000*cos(theta)^4-3835202000*cos(theta)^2+60075480609))/(sin(theta)*cos(theta))}),

but then I get problem with the limit of theta goes to Pi/2.

Third question: is it mathematically legal to say that the vector ['infinity', 0 , 1] is equivalent to the vector [1, 0, 0]? If so, then I think I can become satisfied with my result.

Thank you for the response, Markiyan.

First of all, there should be a zero in the last row, second column of the original matrix I posted. Sorry about that mistake. However, your suggestion solved my problem and allowed me to calculate the eigenvectors. So thank you very much for that tip.

The trigonometric functions still cause some issues in my answer. Here is one of the eigenvectors I get:

u := Vector(3, {(1) = 489794*sin(theta)*cos(theta)/(497206*cos(theta)^2-245103-sqrt(7315644000*cos(theta)^4-3835202000*cos(theta)^2+60075480609)), (2) = 0, (3) = 1})

If I write eval(u[1], theta=0), I get zero as the answer. However, by a quick check I see that both the numerator and denominator are both zero for theta=0. Letting theta go to zero with the help of limit(u[1], theta=0) makes the result 'undefined' (I guess this means answer goes to infinity), as it should be. First question is: why doesn't eval() recognize that there is division by zero?

I have checked by hand that the eigenvectors are correct. Second question: is there some way to avoid this division by zero? I have tried for example to scale the eigenvectors, but that led to a zero-vector. I also tried to "invert" the eigenvector, such that it became

u := Vector(3, {(1) =1, (2) = 0, (3) = (1/489794)*(497206*cos(theta)^2-245103-sqrt(7315644000*cos(theta)^4-3835202000*cos(theta)^2+60075480609))/(sin(theta)*cos(theta))}),

but then I get problem with the limit of theta goes to Pi/2.

Third question: is it mathematically legal to say that the vector ['infinity', 0 , 1] is equivalent to the vector [1, 0, 0]? If so, then I think I can become satisfied with my result.