The expression will involve elliptic integrals, which cannot be expressed in terms of elementary functions. In a computer program, you will most likely have to evaluate it numerically.
Instead, I would recommend combining the approaches: approximate the curve with short line segments.
First, note that the expression for the surface of the oblate spheroid can be rewritten as
$$(R_p x)^2 + (R_p y)^2 + (R_c z)^2 = (R_c R_p)^2$$
If you have some vector $(x, y, z)$, you can use the above to scale it (multiply each of its coordinate) by $c$ to stretch or shrink it to the surface of the oblate spheroid:
$$c = \sqrt{\frac{R_c^2 R_p^2}{(R_p x)^2 + (R_p y)^2 + (R_c z)^2}}$$
This is useful, because you can use the latitude and longitude on an unit sphere to give you the direction, then scale that direction by $c$ to get a point on the oblate spheroid.
Second, if you have point $(x_0, y_0, z_0)$ on the unit sphere ($x_0^2 + y_0^2 + z_0^2 = 1$) corresponding to the latitude and longitude of the starting point on the oblate spheroid, and point $(x_1, y_1, z_1)$ on the unit sphere corresponding to the latitude and longitude of the finishing point on the oblate spheroid, you can calculate the axis of rotation $(x_n, y_n, z_n)$ between these two points using
$$\left\lbrace ~ \begin{aligned}
x_n &= y_0 z_1 - z_0 y_1 \\
y_n &= z_0 x_1 - x_0 z_1 \\
z_n &= x_0 y_1 - y_0 x_1 \\
\end{aligned} \right .$$
which is just a cross product between the two directions. Note that you'll want to normalize this to unit length,
$$\left\lbrace ~ \begin{aligned}
X_n &= \displaystyle \frac{x_n}{x_n^2 + y_n^2 + z_n^2} \\
Y_n &= \displaystyle \frac{y_n}{x_n^2 + y_n^2 + z_n^2} \\
Z_n &= \displaystyle \frac{z_n}{x_n^2 + y_n^2 + z_n^2} \\
\end{aligned} \right.$$
The angular distance between the two directions (assuming unit lengths, $x_0^2 + y_0^2 + z_0^2 = 1$ and $x_1^2 + y_1^2 + z_1^2 = 1$) is $\theta$,
$$\cos\theta = x_0 x_1 + y_0 y_1 + z_0 z_1$$
i.e.
$$\theta = \arccos(x_0 x_1 + y_0 y_1 + z_0 z_1)$$
If you want to split the path into $N$ line segments, you'll need $N+1$ points along the curve on the oblate spheroid between the starting and ending points. These line segments will vary slightly in length, as the oblateness of the spheroid will affect things, but it works very well for approximating the curve length. To get point $k = 0, 1, \dots, N$, use
$$\varphi_k = \frac{k \theta}{N}$$
and apply Rodrigues' rotation formula to the starting point (on unit sphere) to find the $k$'th direction, then scale it to the oblate spheroid. In Cartesian coordinate form, the rotation is
$$\left\lbrace ~ \begin{aligned}
x_k &= x_0 \cos(\varphi_k) + (y_n z_0 - z_n y_0) \sin(\varphi_k) + x_n P \bigr( 1 - \cos(\varphi_k) \bigr) \\
y_k &= y_0 \cos(\varphi_k) + (z_n x_0 - x_n z_0) \sin(\varphi_k) + y_n P \bigr( 1 - \cos(\varphi_k) \bigr) \\
z_k &= z_0 \cos(\varphi_k) + (x_n y_0 - y_n x_0) \sin(\varphi_k) + z_n P \bigr( 1 - \cos(\varphi_k) \bigr) \\
\end{aligned} \right.$$
where $P$ is a constant (dot product between the starting point and the rotation axis),
$$P = x_0 x_n + y_0 y_n + z_0 z_n$$
Note that since $\theta$ is the angular distance, you can make a very sensible choice for $N$ using
$$N = \left\lceil \frac{\theta}{\phi} \right\rceil$$
where $\lceil~\rceil$ denote rounding up (ceil()) and $\phi$ is the largest angular distance you wish to approximate with a straight line. For very small angular distances, you'll get $N = 1$, i.e. just one line segment, or two points, the starting and ending points. Its exact value depends on how precisely you want to calculate the curve length.
Finally, if you have the latitude ($\theta$) and longitude ($\varphi$) of some location on Earth, you can convert those to a point on the unit sphere using
$$\left\lbrace ~ \begin{aligned}
x &= \cos(\theta) \cos(\varphi) \\
y &= \cos(\theta) \sin(\varphi) \\
z &= \sin(\theta) \\
\end{aligned} \right.$$
so that $z$ axis is North (with $\theta \gt 0$ on Northern hemisphere), $x-y$ plane is on the equator, $x$ axis towards zero longitude (zero meridian). If positive $\varphi$ go West and negative East from zero longitude (zero meridian), then positive $y$ axis is in the Americas, and negative $y$ axis in Asia.
It may feel tempting to use $x = R_c \cos(\theta) \cos(\varphi)$ and so on, but that will introduce angular error. Latitude and longitude on Earth are measured as if on an unit sphere. (The error is similar to drawing an ellipse with semimajor axis $2$ and semiminor axis $1$ using $x = 2 \cos\phi$, $y = \sin\phi$. The line through origin to point $(2\cos 45°, \sin 45°)$ does not make a $45°$ angle with the $x$ axis; the angle is only about $26.565°$. Of course, the error in Earth's case is smaller, as the difference in the radiuses is just over $21\text{ km}$, or about $0.34\%$.)
Hope this helps!
