Michael Hardy’s answer gives you a way to find the tangents without derivatives. I’ll expand on N74’s comment as an alternative.
Your ellipse appears to be aligned with the coordinate axes. Wlog, let’s assume that it’s centered on the origin (you can simply translate the resulting lines) so that its equation is $f(x,y)={x^2\over a^2}+{y^2\over b^2}=1$. We have $\nabla f=\left({2x\over a^2},{2y\over b^2}\right)$, so the equation of the tangent line at $(x_0,y_0)$ is, in point-normal form, $$\nabla f(x_0,y_0)\cdot(x-x_0,y-y_0)={2xx_0\over a^2}+{2yy_0\over b^2}-\left({2x_0^2\over a^2}+{2y_0^2\over b^2}\right)={2xx_0\over a^2}+{2yy_0\over b^2}-2=0,$$ or $${xx_0\over a^2}+{yy_0\over b^2}=1.\tag{$*$}$$ Plug the given point into this equation to get a simple linear relationship between $x_0$ and $y_0$, then substitute into the equation of the ellipse to produce a quadratic equation in either $x_0$ or $y_0$.
Another way to find the tangents to any non-degenerate conic that doesn’t involve taking derivatives is to use the dual conic, which consists of all of the lines in the projective plane that are tangent to the conic.
The general equation $Ax^2+Bxy+Cy^2+Dx+Ey+F=0$ of a conic can be written in the form $$\begin{bmatrix}x&y&1\end{bmatrix}\begin{bmatrix}A&\frac B2&\frac D2\\\frac B2&C&\frac E2\\\frac D2&\frac E2&F\end{bmatrix}\begin{bmatrix}x\\y\\1\end{bmatrix}=0$$ or, more compactly, $\mathbf x^TA\mathbf x=0$. If we set $\mathbf x=[x,y,w]^T$, we have the projective equation of the conic. The tangent lines to this conic satisfy the dual equation $\mathscr lA^{-1}\mathscr l^T=0$ (here I’m using row vectors for covectors that represent lines). Since we’re assuming a non-degenerate conic, we know that $A$ is nonsingular. In practice, it’s often more convenient to work with the adjugate matrix $(\det A)A^{-1}$ or some other non-zero multiple of the inverse matrix instead. A point $\mathbf p$ lies on the line $\mathscr l$ if $\mathscr l\mathbf p=0$ (this is the projective equivalent of the point-normal equation of a line). So, finding the tangent lines to a conic that pass through a point $\mathbf p$ amounts to solving the system $$\mathscr lA^{-1}\mathscr l^T=0\\\mathscr l\mathbf p=0.$$ There will, of course be zero, one or two solutions to this system, and each solution will have a free variable, as one might expect when working in homogeneous coordinates.
You appear to be working with ellipses of the form ${(x-h)^2\over a^2}+{(y-k)^2\over b^2}=1$. Expanding this equation and rearranging eventually produces $$A=\begin{bmatrix}b^2&0&-b^2h\\0&a^2&-a^2k\\-b^2h&-a^2k&b^2h^2+a^2k^2-a^2b^2\end{bmatrix}$$ and with a little bit of work we can find a relatively simple expression for the dual conic matrix: $$A^{\tiny\triangle}=\begin{bmatrix}h^2-a^2&hk& h\\hk&k^2-b^2&k\\h&k&1\end{bmatrix}.$$ Letting $\mathbf p=[x_p,y_p,1]^T$ and $\mathscr l=[\lambda,\mu,\nu]$ we get the system of equations $$\mathscr lA^{\tiny\triangle}\mathscr l^T=(h^2-a^2)\lambda^2+(k^2-b^2)\mu^2+\nu^2+2hk\lambda\mu+2h\lambda\nu+2k\mu\nu=0 \\ \mathscr l\mathbf p=x_p\lambda+y_p\mu+\nu=0.$$ The general solution isn’t pretty, but the solutions for specific cases aren’t that hard to work through: since you’re going to be left with a free variable, anyway, you might be able to set one of $\lambda$, $\mu$, or $\nu$ to a convenient value that simplifies the computation.
Then again, if you’re working with the matrix form of the conic equation, you can find the points of tangency by computing the intersection of the polar line of $\mathbf p$ with the conic, i.e., by solving the system of equations $$\mathbf x^TA\mathbf x=0\\\mathbf p^TA\mathbf x=0.$$ The equations for the tangent lines can be found using the usual methods, or by taking advantage of the fact that the polar line of a point on the conic is the tangent at that point.
I’ll sketch out one other method that doesn’t require solving any equations. For the unit circle, it’s very easy to work out geometrically that the points at which lines through $(d,0)$, $d\gt1$ are tangent to it are $\left(\frac1d,\pm{\sqrt{d^2-1}\over d}\right)$. You can scale, rotate and translate the unit circle to produce any ellipse and because affine transformations maintain tangents to circles/ellipses, the same transformation applied to these points will give you the points at which lines through the image of $(d,0)$ are tangent to the ellipse. From there, you can generate the equations of the tangent lines from the pairs of points. Alternatively, you can work in homogeneous coordinates and transform the covectors $\left[-1,\pm\sqrt{d^2-1},d\right]$ of the tangent lines directly, but keep in mind that these are covariant vectors which transform differently than points.
If $P$ is the point through which you’re drawing these tangents, let $Q$ be the inverse image of $P$ under the above transformation. Set $d$ to the distance of $Q$ from the origin, then apply an additional rotation that maps this point to $Q$ before applying the unit-circle-to-ellipse transformation. Given that you’re working with an ellipse in standard position that’s been translated by some amount, this method might work well as the transformation and its inverse are very easy to compute.