Calculating the intersection point of ray vs. arbritarily positioned and rotated ellipsoid - Problem

Question

Let's assume I have an ellipsoid positioned at a given center $\vec{c}$, with principal axes $\vec{e}_1$, $\vec{e}_2$ and $\vec{e}_3$ (calculated by applying three rotation matrices $R_z(\gamma) \cdot R_y(\beta) \cdot R_x(\alpha)$ on each standard basis vector; only the angles are given) and given radii a, b and c. My goal is to calculate the intersection point of a ray (defined as $\vec{p} + t \cdot \vec{q}$) with this ellipsoid (assumed there is one).

My approach so far:

Any ellipsoid can be written as $(\vec{x} - \vec{c})^TA(\vec{x} - \vec{c}) = 1$ with A being a symmetric matrix. According to the principal axis theorem, I can write $A = SDS^T$ with D a diagonal matrix consisting of the eigenvalues of A (in this case $\frac{1}{a^2}, \frac{1}{b^2}, \frac{1}{c^2}$) and S the matrix of the eigenvectors of A (in this case $S = (\vec{e}_1, \vec{e}_2, \vec{e}_3)$). Thus, it follows: $(\vec{x} - \vec{c})^T S\cdot\sqrt{D}\cdot\sqrt{D}\cdot S^T(\vec{x} - \vec{c}) = (\sqrt{D} \cdot S^T \cdot (\vec{x} - \vec{c}))^T (\sqrt{D} \cdot S^T \cdot (\vec{x} - \vec{c})) = 1$. By applying this transformation on my ray, I can define $\vec{p}_{new} = \sqrt{D} \cdot S^T \cdot (\vec{p} - \vec{c})$ and $\vec{q}_{new} = \sqrt{D} \cdot S^T \cdot \vec{q}$ such that I can calculate the intersection of a new ray $\vec{p}_{new} + t \cdot \vec{q}_{new}$ with a unit sphere around the origin. This is how I implemented it in Matlab. The resulting variable t can then be used for finding the intersection point with the initial ray.

Problem: There seems to be nothing wrong with the approach itself, as for y- and z-rotations this works fine. However, it doesn't calculate x-rotations properly (to be more specific, it doesn't apply any x-rotations at all). I assume there is something wrong with the way I order my eigenvalues and eigenvectors but I'm not really sure about this so any tip on what I might be missing here is greatly appreciated.

Code:

        Rx=rotx(alpha); Ry=roty(beta); Rz=rotz(gamma);
        x_new = Rz*Ry*Rx*[1 0 0]';
        y_new = Rz*Ry*Rx*[0 1 0]';
        z_new = Rz*Ry*Rx*[0 0 1]';
        D = [1/a 0 0; 0 1/b 0; 0 0 1/c];
        S = [x_new(1) y_new(1) z_new(1); x_new(2) y_new(2) z_new(2); x_new(3) y_new(3) z_new(3)];
        M = D*transpose(S);
        p_new = M*([p(1) p(2) p(3)]' - [c(1) c(2) c(3)]');
        q_new = M*([q(1) q(2) q(3)]');
        F(1)=q_new(1)^2+q_new(2)^2+q_new(3)^2;
        F(2)=2*(p_new(1))*q_new(1)+2*(p_new(2))*q_new(2)+2*(p_new(3))*q_new(3);
        F(3)=(p_new(1))^2+(p_new(2))^2+(p_new(3))^2-1;
        s1=roots(F);