How do you find a geodesic on an ellipsoid?
So a geodesic is essentially the minimization of a functional, that is, a function of functions. They are the main study of Calculus of Variations. Now, if we have a functional of the form $$J[y]=\int_a^b F(x,y,y') \ dx$$ defined on the set of functions $y(x)$ which have continuous first derivatives in $[a,b]$, and satisfy the boundary conditions $y(a)=A$, $y(b)=B$, then a neccesary condition for $J[y]$ to have an extremum for a function $y(x)$ is that $y(x)$ satisfy the Euler-Lagrange equation $$F_y-\frac{d}{dx}F_{y'}=0.$$
In essence, we need to solve that above differential equation.
Without getting anymore nitty gritty on things, the functional we care about here is the arc length integral.
Firstly, suppose we have a surface $\mathcal{S}$ defined by the vector equation $$\vec{r}(u,v)=x(u,v) \hat{i} + y(u,v) \hat{j} + z(u,v) \hat{k}.$$ The geodesic curve lying on surface $\mathcal{S}$ can be specified by the equations $$u=u(t), \qquad v=v(t),$$ and can be found by minimizing the arc length integeral $$L=\int_{t_0}^{t_1} \sqrt{\left( \frac{dx}{dt} \right)^2+\left( \frac{dy}{dt} \right)^2+\left( \frac{dz}{dt} \right)^2} \ dt$$ But since $x,y,z$ are each functions of more than one variable, the chain rule is going to turn this integral into a monstrosity. Thankfully, it's a monstrosity with a pretty neat pattern so we can make things look a bit prettier.
Letting $u'=\frac{du}{dt}$ and $v'=\frac{dv}{dt}$, the integeral above can just be rewritten as
$$J[u,v]=\int_{t_0}^{t_1} \sqrt{Eu'^2+2Fu'v'+Gv'^2} \ dt,$$
where $E$, $F$, and $G$ are the coefficients of the first fundamental (quadradic) form of the surface, i.e., $$\begin{cases} \displaystyle{E=\vec{r}_u\cdot \vec{r}_u = \left( \frac{\partial x}{\partial u} \right) ^2 +\left( \frac{\partial y}{\partial u} \right) ^2 + \left( \frac{\partial z}{\partial u} \right) ^2} \\ \\ \displaystyle{F=\vec{r}_u \cdot \vec{r}_v = \frac{\partial x}{\partial u} \frac{\partial x}{\partial v} + \frac{\partial y}{\partial u} \frac{\partial y}{\partial v} + \frac{\partial z}{\partial u} \frac{\partial z}{\partial v}}\\ \\ \displaystyle{G=\vec{r}_v\cdot \vec{r}_v=\left( \frac{\partial x}{\partial v} \right) ^2 +\left( \frac{\partial y}{\partial v} \right) ^2 + \left( \frac{\partial z}{\partial v} \right) ^2}. \end{cases}$$
The Euler-Lagrange equation in this case corresponds to the two different equations $$\displaystyle{ F_u- \frac{d}{dt} F_{u'} = 0,} \ \ \ \ \ \ \ \displaystyle{F_v- \frac{d}{dt} F_{v'} = 0,}$$ hence, we obtain \begin{equation}\label{geodesiceq1} \frac{E_u u'^2+2F_u u' v' + G_u v'^2}{\sqrt{Eu'^2+2Fu'v'+Gv'^2}} - \frac{d}{dt} \frac{2(Eu'+Fv')}{\sqrt{Eu'^2+2Fu'v'+Gv'^2}} = 0, \end{equation} \begin{equation}\label{geodesiceq2} \frac{E_v u'^2+2F_v u' v' + G_v v'^2}{\sqrt{Eu'^2+2Fu'v'+Gv'^2}} - \frac{d}{dt} \frac{2(Fu'+Gv')}{\sqrt{Eu'^2+2Fu'v'+Gv'^2}} = 0. \end{equation}
Which are the two differential equations whose solutions provide the geodesic on surface $\mathcal{S}$.
Now let me note one important case. The case where $E$ and $G$ are explicit functions of $v$ only and $F=0$, we have $$\frac{E_v+v'^2G_v}{2\sqrt{E+Gv'^2}}-\frac{d}{du} \left( \frac{Gv'}{\sqrt{E+Gv'^2}} \right) =0,$$ so $$ E_v + v'^2G_v - 2 G \sqrt{E+Gv'^2} \bigg[ \frac{v''}{\sqrt{E+Gv'^2}} + \left( \frac{1}{2} \right) \frac{v'(2Gv'v'')}{(E+Gv'^2)^{\frac{3}{2}}} \bigg] =0 $$ $$E_v + v'^2 G_v - 2Gv'' + \frac{2G^2v'^2v''}{E+Gv'^2} = 0 $$ \begin{equation} \frac{Gv'^2}{\sqrt{E+Gv'^2}} - \sqrt{E+Gv'^2} = c_1 \end{equation} which can be made even more helpful by noting $\displaystyle{v'= \frac{dv}{du}}$, giving us $$Gv'^2-(E+Gv'^2)=c_1 \sqrt{E+Gv'^2}$$ $$\left( - \frac{E}{c_1} \right) ^2 = E + Gv'^2$$ $$\frac{E^2-{c_1}^2E}{G{c_1}^2} = v'^2,$$ finally providing \begin{equation} u={c_1} \int_{v_1=v(u_1)}^{v_2=v(u_2)} \sqrt{ \frac{G}{E^2-{c_1}^2E }} \ dv. \end{equation}
Now lets talk about the ellipsoid. An ellipsoid with semi-axis lengths $a,b,c$ is given by the vector equation $$\vec{r}(u,v)= \langle a \cos(u) \sin(v), b \sin(u) \sin(v), c \cos(v) \rangle$$ with $u \in [0,2\pi)$ and $v\in [0, \pi]$.
The first fundamental form of the ellipsoid which I sadly found by hand before I realized I could've went to Google for is
$$\begin{cases} E= b^2 \cos^2(u) \sin^2(v)+a^2 \sin^2(u) \sin^2(v) \\ F= \left( b^2-a^2 \right) \sin(v) \cos(v) \sin(u) \cos(u)\\ G= \cos^2(v) \left[ (a^2-b^2) \cos^2(u)+b^2 \right] +c^2 \sin^2(v) \end{cases}$$
Now listen, I'm actually procrasinating the hell out of a transfer essay by being here and looking at the fundamental form above is pretty worrying, because it's far from elegant and being easy to deal with and would entail solving that frightening differential equation I mentioned earlier.
So I'm going to let $a=b$. Playing around with $c$ is ellipsoid-y enough for me, and if you want to go for the full $a,b,c$ then be my guest, haha! So doing what I said, things are a lot prettier:
$$\begin{cases} E= a^2 \sin^2(v) \\ F= 0 \\ G=b^2 \cos^2(v) + c^2 \sin^2(v) \end{cases}$$
and this is where we can apply that last remark I made earlier regarding the first fundamental form looking like this. We thus have the geodesic on this surface given by
$$u={c_1} \int \sqrt{ \frac{G}{E^2-{c_1}^2E }} \ dv = c_1 \int \sqrt{\frac{b^2 \cos^2(v) + c^2 \sin^2(v)}{a^4 \sin^4(v) - c_1^2 a^2 \sin^2(v)}} \ dv.$$
Now in the case of a sphere, solving the integeral and rearranging gives a plane -- and that planes interesection with the sphere represents geodesics on a sphere. I myself couldn't get anywhere with integeral, but I'm sure someone else can so I'll just leave things off here.
... Now back to my essay.
On a very final note, here is my research paper from my first semester. I prove every step up to proving the very Euler-Lagrange equations themselves and then show some examples of it in use, like finding the geodesic on a sphere and geodesic of a function revolved around an axis.
Like Jack mentioned in the comments to your post, if you consider just the ellipsoids that can be obtained through revolving a function $f(x)$, the integeral is much easier. I prove it in the above paper, but here it is again:
$$v=c_1 \int_{u_1}^{u_2} \frac{\sqrt{1+(f'(u))^2}}{f(u) \sqrt{(f(u))^2-{c_1}^2}} \ du$$
Note, $x=u$ in the above integeral.
The following addresses only initial values only.
If all geodesics are tangent to circle
$$ x^2+y^2= r_{min}^2 @ \,r=r_{min}$$
as the Clairaut constant then for an ellipsoid of revolution
$$ (z/c)^2+(r/a)^2 = (z/c)^2+(x^2+y^2)/a ^2 = 1 $$
we differentiate to find meridian slope
$$ \tan \phi=-(a/c) \sqrt{(a/r)^2-1}, $$
and $r=f(\theta) $ is found by integration (to get elliptic integral result)
$$ \frac{dr}{r \, d\theta}= \sin\phi \, \sqrt{1-(r/r_{min})^2 }.$$
similarly $z$ can be set up. We have $ (r,z)$ as functions of $\theta$
With given $ (r_B,r_A), \theta_{AB}, (z_B,z_A) $ we cannot directly use this, however..
we can start with initial conditions $(r_B, \theta_A =0, z_A ) $ but to arrive at last point we can vary /adjust Clairaut constant $r_{min}$ several small increments in a numerical procedure.