Points of tangency from a point M to an ellipsoid-plane intersection Having an spheroid $S$ of semi-major axis $a$ in the equatorial plane in direction $x-axis$ and $y-axis$, and of semi-minor axis $b$ in direction $z-axis$.
$(S): (\frac{x}{a})^2+(\frac{y}{a})^2+(\frac{z}{b})^2=1$
A plane $(P): ux+vy+wz=d$ intersects $S$
That yields an ellipse $E$ in 3D.
I Intersected $S$ and $P$, I got a quadratic equation of quadratic conics, in the form:
$(E): Ax^2+Bxy+Cy^2+Dx+Ey+F=0$
After equating and solving, I got the coefficients of the ellipse $E$
$A=(wb)^2+(ua)^2$
$B=2uv(a)^2$
$C=(wb)^2+(va)^2$
$D=-2du(a)^2$
$E=-2dv(a)^2$
$F=a^2(d^2-(bw)^2)$
My question is as follow:
I have a point $M(i,j,k)$ lying in the plane of intersection, and it is outside the ellipse $E$, how to find the coordinates of the two tangent points $B$ and $C$ from the point $M$ to the ellipse $E$?
Please see figure from geogebra:

Add

 A: Equation of the ellipsoid:
$$\frac{x^2+y^2}{a^2}+\frac{z^2}{b^2}=1 \tag{1}$$
Equation of the polar plane w.r.t. $M$:
$$\frac{ix+jy}{a^2}+\frac{kz}{b^2}=1 \tag{2}$$
Equation of the plane $P$ containing $M$:
$$ux+vy+wz=d \tag{3}$$
under the constraint $ui+vj+wk=d$.
In principle, $(1) \cap (2) \cap (3)$, which is the required intersections, is solvable with discriminant
$$\frac{i^2+j^2}{a^2}+\frac{k^2}{b^2}-1>0$$
A: First of all, the equation that you’ve derived ostensibly by solving equation (P) for $z$ and then substituting into (S) isn’t the equation of the intersection curve. Indeed, a curve in a 3-D space can’t be represented by a single implicit Cartesian equation. You’ve got the equation of an elliptical cylinder that happens to have the same intersection with the cutting plane. The intersection of this cylinder with the $x$-$y$ plane is the orthogonal projection of the ellipse, so is in most cases foreshortened.  
There are several ways to solve this problem, and the best/easiest really depends on what else you want to know about this configuration. If all you want are the two points of tangency, then a straightforward way to go is to find the intersection of $P$ and the polar plane to $M$, then computing the intersection of this line with the spheroid. If you use a standard parametric representation for the line, this is a straightforward matter of solving a single-variable quadratic equation.  
Specifically, the polar plane to $M$ is ${ix\over a^2}+{jy\over a^2}+{kz\over b^2}=1$, so its intersection with $P$ has direction vector $$\mathbf v = \left(\frac i{a^2},\frac j{a^2},\frac k{b^2}\right)\times(u,v,w)$$ and with a bit of work you can find a point on this line, such as $$\mathbf p = \left({dj-a^2v \over ju-iv} , {di-a^2u \over iv-ju}, 0\right).$$ Plug $\mathbf p+t\mathbf v$ into the equation of the spheriod and solve for $t$.  
A related approach is to transform the problem into a simpler one that involves intersection with and tangents to the unit sphere, solve that problem, then transform back. The required transformations are simple scaling operations that preserve intersections and tangency. Mapping the spheroid to the unit sphere takes $(i,j,k)$ to $(i/a,j/a,k/b)$ and takes $P$ to the plane $P'$ with equation $aux+avy+bwz=d$. The points of tangency can be found using the method described above, or you could further rotate or reflect to make $P'$ parallel to the $x$-$y$ plane and then solve the simple problem of finding the ends of the chord of contact to a circle. This additional rotation/reflection also makes it fairly easy to compute a parameterization of the intersection curve: it’s just the unit circle rotated/reflected and scaled back into position. Once you have this parameterization of the intersection curve, you can use it to generate test points for the search algorithm that you’ve described.  
Finally, there’s a version of the tried-and-true method that you likely learned in planar analytic geometry. Generate a one-parameter family of lines through $M$ that lie in $P$ and work out their intersection with the spheriod. For tangency, you want the resulting equation to have a double root, which will then generate a simple quadratic equation in the line parameter. Once you have that, finding the points of tangency will be a matter of back-substitution.
