Traduce el elipsoide al origen, restando $(x_3, y_3, z_3)$ de $P1$ y $P2$, por lo que la línea está parametrizada como $$\vec{p} = \vec{v}_0 + t \vec{v}_1$$ donde $$\begin{aligned} \vec{v}_0 = (\chi_0, \gamma_0, \zeta_0) &= (x_1 - x_3, y_1 - y_3, z_1 - z_3) \\ \vec{v}_1 = (\chi_1, \gamma_1, \zeta_1) &= (x_2 - x_1, y_2 - y_1, z_2 - z_1) \\ \end{aligned}$$ y el elipsoide es $$\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{b^2} = 1$$ Sustituyendo $\vec{p}$ en el elipsoide se obtiene $$\left(\frac{\chi_0 + t \chi_1}{a}\right)^2 + \left(\frac{\gamma_0 + t \gamma_1}{b}\right)^2 + \left(\frac{\zeta_0 + t \zeta_1}{c}\right)^2 = 1$$ Expande los términos y recoge los coeficientes para las potencias de $t$ y obtendrás $$\begin{aligned} \left( \displaystyle \frac{\chi_1^2}{a^2} + \frac{\gamma_1^2}{b^2} + \frac{\zeta_1^2}{c^2} \right) & t^2 ~ + \\ \left( \displaystyle \frac{2 \chi_0 \chi_1}{a^2} + \frac{2 \gamma_0 \gamma_1}{b^2} + \frac{2 \zeta_0 \zeta_1}{c^2} \right) & t ~ + \\ \left( \frac{\chi_0^2}{a^2} + \frac{\gamma_0^2}{b^2} + \frac{\zeta_0^2}{c^2} - 1 \right) & ~ = 0 \\ \end{aligned}$$ Esta es una ecuación cuadrática simple en t, que puede tener 0, 1 o 2 raíces reales $t$.
Recuerda que solo hemos traducido el sistema de coordenadas para que fuera relativo al centro del elipsoide, pero no escalamos ni rotamos el sistema de coordenadas. Por lo tanto, después de haber encontrado $t$, puedes encontrar el punto al que se refiere, en el sistema de coordenadas original, usando tu línea parametrizada original, $$P = P1 + t ( P2 - P1 )$$ es decir $$\left\lbrace ~ \begin{aligned} x &= x_1 + t \chi_1 \\ y &= y_1 + t \gamma_1 \\ z &= z_1 + t \zeta_1 \\ \end{aligned} \right.$$ para cada punto de intersección $t$.