Hamed’s answer covers your specific question in great detail, but I’d like to suggest a couple of different approaches to the problem you were trying to solve that led to it that might make the algebra easier.
Every ellipse can be obtained from the unit circle via an affine transformation, which preserves its tangents. So, start by solving the simpler problem of finding tangents to the unit circle. For the point $P(d,0)$, $d\ge1$, this is particularly simple. Its polar line is $x=1/d$, which intersects the unit circle at $y=\pm\sqrt{1-1/d^2}$. For a general point $(\xi',\eta')$ that is not interior to the circle, set $d=\sqrt{\xi'^2+\eta'^2}$, the distance of the point from the origin, and rotate the two intersection points derived previously for the simple case into place, producing the points $$\frac1{d^2}
\left(\xi'\pm\eta'\sqrt{d^2-1},\eta'\mp\xi'\sqrt{d^2-1}\right).\tag1$$ Using the point/normal form of equation for a line and expanding $d^2$, these points give us the equations $$\left(\xi'\pm\eta'\sqrt{\xi'^2+\eta'^2-1}\right)\,x'+\left(\eta'\mp\xi'\sqrt{\xi'^2+\eta'^2-1}\right)\,y'=\xi'^2+\eta'^2\tag2$$ for the tangent lines to the unit circle.
For an ellipse in standard position with semi-axis lengths $a$ and $b$ and point $P(\xi,\eta)$, transform into the simple case, solve it, and transform back. The transformation to the unit circle is a simple scaling: $\xi'=\xi/a$ and $\eta'=\eta/b$. The transformation back can be accomplished by replacing $x'$ by $x/a$ and $y'$ by $y/b$. Making these substitutions into equation (2) rearranging and multiplying through by $ab$ results in the equations $$b^2\xi(x-\xi)+a^2\eta(y-\eta)\pm(\eta x-\xi y)\sqrt{b^2\xi^2+a^2\eta^2-a^2b^2}=0\tag3$$ for an ellipse in standard position.
Another way to go is to use the dual conic to the ellipse. First, translate the origin to $P(\xi,\eta)$. The matrix of the ellipse in this coordinate system is $$Q=\begin{bmatrix}\frac1{a^2} & 0 & \frac\xi{a^2} \\ 0 & \frac1{b^2} & \frac\eta{b^2} \\ \frac\xi{a^2} & \frac\eta{b^2} & {\xi^2\over a^2}+{\eta^2\over b^2}-1\end{bmatrix}$$ and the matrix of its dual is $$Q^{-1}=\begin{bmatrix}a^2-\xi^2 & -\xi\eta & \xi \\ -\xi\eta & b^2-\eta^2 & \eta \\ \xi & \eta & -1 \end{bmatrix}.$$ The tangents we’re interested in pass through the origin, so will have the homogeneous representation $[\lambda,\mu,0]$, which is equivalent to zeroing out the last row and column of the above dual matrix. To get the equation of the two lines, we invert the upper-right $2\times2$ submatrix. Since we’re working in homogeneous coordinates, any nonzero scalar multiple of this matrix is equivalent, so we use the adjugate instead to make this computation trivial, resulting in the matrix $$T=\begin{bmatrix} b^2-\eta^2 & \xi\eta & 0 \\ \xi\eta & a^2-\xi^2 & 0 \\ 0&0&0 \end{bmatrix}$$ for the degenerate conic that consists of the tangents to the ellipse through the origin. Finally, translate $P$ back to where it started, yielding as an equation of the tangent lines through $P$ $$(b^2-\eta^2)(x-\xi)^2 + 2\xi\eta(x-\xi)(y-\eta) + (a^2-\xi^2)(y-\eta)^2 = 0$$ which simplifies to $$b^2(x-\xi)^2+a^2(y-\eta)^2-(\eta x-\xi y)^2=0.\tag4$$ (If you multiply the two equations in (3) together and factor, you’ll get equation (4) times a positive constant.)
If you need the individual equations, you can painstakingly factor equation (4), cheat and use equation (3), or go back to the dual and solve $[\lambda,\mu,0]Q^{-1}[\lambda,\mu,0]^T=0$ for $\lambda$ and $\mu$: $$(a^2-\xi^2)\lambda^2 - 2\xi\eta\lambda\mu + (b^2-\eta^2)\mu^2 = 0$$ which yields as one possible solution the lines $$\left[b^2-\eta^2,\xi\eta\pm\sqrt{b^2\xi^2+a^2\eta^2-a^2b^2},0\right]\tag5$$ in homogeneous form, or $$(b^2-\eta^2)\,x+\left(\xi\eta\pm\sqrt{b^2\xi^2+a^2\eta^2-a^2b^2}\right)\,y=0$$ in Cartesian, which after translation becomes $$(b^2-\eta^2)(x-\xi)+\left(\xi\eta\pm\sqrt{b^2\xi^2+a^2\eta^2-a^2b^2}\right)(y-\eta)=0.\tag6$$ Unfortunately, this last equation gives only one of the tangent lines when $\eta=\pm b$, which will have to be handled as a fairly straightforward special case that I’ll leave to you. Equation (4), on the other hand, is correct in all cases. This method generalizes to higher dimensions, allowing the tangent cone to an ellipsoid to be found the same way, for instance.
Additional method: In this answer, MvG describes a method to find the intersection of a line and conic that involves decomposing a rank 2 matrix into its two generating vectors. Since this intersection problem is the dual of the original problem you were trying to solve, that of finding the tangents to an ellipse through a given point, the same method can be used to find the solution to that original problem directly.
This time we’ll leave the origin at the center of the ellipse. Its matrix is then simply $C=\operatorname{diag}(1/a^2,1/b^2,-1)$. The matrix $P_\times^TC^{-1}P_\times$, where $P_\times$ is the “cross-product matrix” for the point $P$, is a degenerate conic that consists of the two tangent lines through $P$: If the point $X$ lies on this conic, then $$X^TP_\times^TC^{-1}P_\times X=(P\times X)^TC^{-1}(P\times X)=0$$ which means that the line $P\times X$ through $X$ and $P$ is tangent to the ellipse.
For $P(\xi,\eta)$ we have $$P_\times=\begin{bmatrix}0&-1&\eta\\1&0&-\xi\\-\eta&\xi&0\end{bmatrix}$$ and $C^{-1}$ is trivial to compute, so for the tangent lines we end up with the matrix $$M=\begin{bmatrix}b^2-\eta^2 & \xi\eta & -b^2\xi \\\xi\eta & a^2-\xi^2 & -a^2\eta \\ -b^2\xi & -a^2\eta & a^2\eta^2+b^2\xi^2 \end{bmatrix}.$$ (Observe that this is the matrix $T$ above translated by $-(b^2\xi,a^2\eta)$. Observe also that $T=O_\times^TQ^{-1}O_\times$, where $O$ is the origin. If you examine $O_\times$, you’ll see why we could use the adjugate of the upper-left submatrix in that computation.) It’s fairly easy to verify that the null vector of this matrix is the point of intersection, $P(\xi,\eta)$.
At this point we can say that we’re done since we’ve got an equation for the two lines—$\mathbf x^TM\mathbf x=0$—but we’ll continue and decompose $M$ into its constituent lines $\mathbf l$ and $\mathbf m$. $M$ is a multiple of $\mathbf l\mathbf m^T+\mathbf m\mathbf l^T$, so following the method described in the linked answer, we form $M+\lambda P_\times$, set some $2\times2$ minor to zero and solve for $\lambda$. The upper-left minor looks like the simplest to work with, so we have $$(a^2-\xi^2)(b^2-\eta^2)-(\xi^2\eta^2-\lambda^2)=0$$ from which $$\lambda=\pm\sqrt{b^2\xi^2+a^2\eta^2-a^2b^2}.$$ We’ll use the positive root and then take for $\mathbf l$ and $\mathbf m$ any convenient row and column of $M+\lambda P_\times$, yielding the two lines $$\left[b^2\xi\pm\eta\sqrt{b^2\xi^2+a^2\eta^2-a^2b^2},a^2\eta\mp\xi\sqrt{b^2\xi^2+a^2\eta^2-a^2b^2},-(b^2\xi^2+a^2\eta^2)\right].$$ Converting these vectors into Cartesian equations and rearranging a bit produces equation (3).
Observe that this last method is completely general. Having the ellipse in standard position made the calculations simpler, but the result doesn’t depend on this. Neither does the method depend on the conic’s being an ellipse. It will work just as well with a parabola or hyperbola.
