The same maximal distance occurs in each quadrant, so we can restrict attention to $t\in[0,\pi/2]$. The tangent vector at $t$ is $(-a\sin t,b\cos t)$. This vector is normal to the line, so we just have to take the scalar product of a unit vector in this direction with the position vector in order to find the distance of the origin from the line:
$$
\begin{eqnarray}
D
&=&
\left\lvert\frac{(-a\sin t,b\cos t)}{\sqrt{a^2\sin^2t+b^2\cos^2t}}\cdot(a\cos t,b\sin t)\right\rvert
\\
&=&
\frac{(a^2-b^2)\sin t\cos t}{\sqrt{a^2\sin^2t+b^2\cos^2t}}\;.
\end{eqnarray}$$
Differentiating with respect to $t$ yields
$$\frac{a^2\sin^4 t-b^2\cos^4t}{\left(a^2\sin^2t+b^2\cos^2t\right)^{3/2}}\;,$$
and setting this to zero yields
$$a^2\sin^4t=b^2\cos^4t\;,$$
$$t=\arctan\sqrt{\frac ba}\;.$$
Using $\cos t=1/\sqrt{1+\tan^2 t}$, we can evaluate $D$ at this parameter:
$$
\begin{eqnarray}
D
&=&
\frac{(a^2-b^2)\sin t\cos t}{\sqrt{a^2\sin^2t+b^2\cos^2t}}
\\
&=&
\frac{(a^2-b^2)\tan t}{\sqrt{a^2\tan^2t+b^2}}\cos t
\\
&=&
\frac{(a^2-b^2)\tan t}{\sqrt{a^2\tan^2t+b^2}}\frac1{\sqrt{1+\tan^2 t}}
\\
&=&
\frac{(a^2-b^2)\sqrt{b/a}}{\sqrt{a^2(b/a)+b^2}}\frac1{\sqrt{1+b/a}}
\\
&=&
\frac{a^2-b^2}{a+b}\;.
\\
&=&
a-b\;.
\end{eqnarray}$$
The result obviously supports your idea that there might be a simpler way to do this.
Which of the two $k$ values do I use in solving for $r'$ and $z'$ ?
Either one will be a point of intersection, since in general a line intersects an ellipse in two points. Which one you need, or if it doesn't matter, depends on your application.
How can you tell when the line doesn't intersect the ellipse? Is the quadratic equation for $k$ not have any real roots?
Right. And if you only have a single root, i.e. both roots coincide, then the line will be a tangent.
Is it important for $r$ to be $>0$ in for this equation even though I'm not testing a point, but a line instead?
I see no reason to require this here.
Why is the semi-minor axis being defined as $a_e(1-f)$? I usually just define the semi-minor axis the same way I define the major, so could I just replace all the $(1-f)$s in the above equations with my desired semi-minor axis length?
Yes, you can write the semi-minor axis (called $a_p$ in your document) instead of $a_e(1-f)$. Which means you could replace the $(1-f)$ themselves by $\frac{a_p}{a_e}$ or multiply the whole equation by $a_e^2$.
Finally, is there any simpler, faster way to see if and were an ellipse and line intersect?
Depends on how your objects are given. If you really have an ellipse and a line in the described form, I can't think of something fundamentally easier. If, on the other hand, you have an ellipse in some other representation, then you might be able to use that representation directly to compute the intersection, instead of transforming its representation first. In any case, you can't avoid the quadratic equation.
Best Answer
I'd follow this answer of mine.
Turn your conics into symmetric matrices: $$A=\begin{bmatrix}4&-1&0\\-1&2&1\\0&1&-8\end{bmatrix}\qquad B=\begin{bmatrix}2&0&0\\0&-1&0\\0&0&-1\end{bmatrix}$$ Note that I scaled your first equation by a factor of $2$ because symmetrically distributing the off-diagonal terms divides them by two and I wanted to avoid fractions.
Find linear combinations with zero determinant. $$\det(\lambda A+\mu B)=-60\lambda^3-9\lambda^2\mu+16\lambda\mu^2+2\mu^3$$ This does not factor so you can't expect rational solutions. In particular zero is not a solution for either $\lambda$ or $\mu$ (although $\lambda=\mu=0$ is a solution but we ignore that). So set $\lambda=1$ and compute $\mu$ for that: $$\mu_1\approx-8.10\qquad\mu_2\approx-1.88\qquad\mu_3\approx1.97$$ The exact solutions are tedius to write out.
$A+\mu_iB$ is a degenerate conic formed by a pair of lines passing through the points of intersection. I'll continue with $\mu_1$ for my numeric examples. $$C_1=A+\mu_1B\approx\begin{bmatrix}-12.2&-1&0\\-1&10.1&1\\0&1&0.0982\end{bmatrix}$$
The adjugate of that is $$\operatorname{adj}C_1\approx\begin{bmatrix}-0.00805&0.0982&-1\\0.0982&-1.20&12.2\\-1&12.2&-124\end{bmatrix}$$ It has rank 1, so all rows and all columns are multiple of one another. Take any non-zero row or column of that, e.g. the first row, and you have the homogeneous coordinates where the two lines intersect. (Divide by the last coordinate if you want to have regulas $(x,y)$ coordinates but you don't need that.) Now use these coordinates to form an anti-symmetric matrix $$P_1\approx\begin{bmatrix}0&-1&-0.0982\\1&0&-0.00805\\0.0982&0.00805&0\end{bmatrix}$$ Then consider $$C_1+\lambda P_1\approx\begin{bmatrix}-12.2&-1-\lambda&-0.0982\lambda\\-1+\lambda&10.1&1-0.00805\lambda\\0.0982\lambda&1+0.00805\lambda&0.0982\end{bmatrix}$$ Take any $2\times2$ minor of this, e.g. $$\begin{vmatrix}-12.2&-1-\lambda\\-1+\lambda&10.1\end{vmatrix}\approx\lambda^2-124\overset!=0$$ From this you can conclude that $\lambda\approx\pm11.1$. Choosing the positive solution (an arbitrary choice) you get $$C_1+\lambda P_1\approx\begin{bmatrix}-12.2&-12.1&-1.09\\10.1&10.1&0.910\\1.09&1.09&0.0982\end{bmatrix}$$ This matrix has again rank one, so its rows are multiples of one another, and its columns are multiples of one another. Pick any non-zero row and any non-zero column, and you have the equations of two of the lines. $$g_1\approx[-12.2:-12.1:-1.09]\qquad h_1\approx[-12.2:10.1:1.09]$$.
Repeat these steps for $\mu_2$. I wont print this here, but in the end you obtain $$g_2\approx[0.248:-1.20:1.23]\qquad h_2\approx[0.248:-0.799:-1.23]$$ Now intersect one of the lines from roun $1$ with one from round $2$ by computing the corss product to obtain one of your points of intersection. Divide by the last coordinate to dehomogenize. \begin{align*}q_1&=g_1\times g_2\approx[-16.3:14.7:17.7]\to(-0.921, 0.835)\\q_2&=g_1\times h_2\approx[14.1:-15.3:12.8]\to(1.10, -1.20)\\q_3&=h_1\times g_2\approx[13.8:15.3:12.1]\to(1.14, 1.26)\\q_4&=h_1\times h_2\approx[-11.6, -14.7, 7.23]\to(-1.61, -2.04)\end{align*}