It is a bit more difficult than the formula for circles.
You can do the following (which can also easily be implemented into an algorithm):
- Find the center points $c_1,c_2$ of each ellipse
- Define the line segment between those center points
- Find the intersection of this segment with each ellipse, denoted by $p_1$ and $p_2$
- The distance you're looking for is $|p_1-p_2|$
This can be solved in closed form, here the important steps to get the algorithm running:
Center of ellipses are $c_1 = (x_1,y_1)$ and $c_2=(x_2,y_2)$, therefore the line segment between those points is given as
\begin{align}
y & = \frac{y_2-y_1}{x_2-x_1} (x-x_1) + y_1 \\
&=\frac{y_2-y_1}{x_2-x_1} (x-x_2) + y_2
\end{align}
Let us now calculate the intersection with the first ellipse $E_1: \frac{(x-x_1)^2}{a_1^2} +\frac{(y-y_1)^2}{b_1^2}=1$ by just inserting the formula for $y$:
\begin{align}
\frac{(x-x_1)^2}{a_1^2} +\frac{(y-y_1)^2}{b_1^2}&=1 \\
\frac{(x-x_1)^2}{a_1^2} +\frac{(\frac{y_2-y_1}{x_2-x_1} (x-x_1) + y_1-y_1)^2}{b_1^2}&= 1 \\
\frac{(x-x_1)^2}{a_1^2} +\frac{(y_2-y_1)^2}{(x_2-x_1)^2b_1^2}(x-x_1)^2&= 1 \\
\left(\frac{1}{a_1^2} +\frac{(y_2-y_1)^2}{(x_2-x_1)^2b_1^2}\right)(x-x_1)^2&= 1 \\
x^{(1)}_{1,2} = \pm \left(\frac{1}{a_1^2} +\frac{(y_2-y_1)^2}{(x_2-x_1)^2b_1^2}\right)^{-\frac 12}+x_1
\end{align}
Remark: $x^{(1)}_{1/2}$ represents the two solutions (subscript $1,2$) of the intersection of the first ellipse with the connection line.
Analogously by entering the other formula for $y$ in the second equation we get
$$
x^{(2)}_{1,2} = \pm \left(\frac{1}{a_2^2} +\frac{(y_2-y_1)^2}{(x_2-x_1)^2b_2^2}\right)^{-\frac 12}+x_2
$$
Because the line intersects with each ellipse at two points we have two different solutions. To find out the correct solution we can make a case-by-case analysis:
If $x_1 \leq x_2$: The $x$-value of the intersection must be inside the interval $[x_1,x_2]$, therefore choose $x^{(1)}$ and $x^{(2)}$ s.t. $x_1 \leq x^{(1)} \leq x^{(2)} \leq x_2$.
If $x_1 > x_2$: Change the order and choose those solutions s.t. $x_1 > x^{(1)} > x^{(2)} > x_2$.
The corresponding $y$ values can be calculated using the formula on top of this answer.
One could calculate the difference explicitly, but that requires more computing and simplifying terms. As you only need an algorithmic approach I think that this is the way to go.