How does one prove that on $S^n$ (with the standard connection) any geodesic between two fixed points is part of a great circle?

For the special case of $S^2$ I tried an naive approach of just writing down the geodesic equations (by writing the Euler-Lagrange equations of the length function) and solving them to gain some insights but even if the equations are solvable I can't see how to show that they are great circles. (the solutions are some pretty complicated functions which don't give me much insight)

I checked the article on Great Circles on Wolfram Mathworld for a coordinate geometry approach to it but that article looked quite cryptic to me!

One knows that on compact semi-simple lie groups any one-parameter subgroup generates a geodesic and $S^n$ is the quotient of 2 compact semi-simple lie groups $SO(n+1)/SO(n)$. Is this line of thought useful for this question?

================================================================================= After some of the responses came let me put in "a" way of seeing the above for $S^2$ (wonder if it is correct). If $\theta$ and $\phi$ are the standard coordinates on $S^2$ then the equations for the curve are

$$\ddot{\theta} = \dot{\phi}^2 sin(\theta)cos(\theta)$$ $$\dot{\phi}sin^2{\theta} = k$$

where $k$ is some constant set by the initial data of the curve.

Now given the initial point I can choose my coordinate system such that the the initial data looks like $\dot{\phi}=0$, $\theta = \text{some constant}$, $\dot{\theta}=\text{some constant}$, $\phi = \text{some constant}$. Then the differential equations tell me that the $k=0$ and the only way it can happen for times is by having ,

$$\dot{\phi} = 0$$

Which clearly gives me a longitude in this coordinate system. Hence the geodesic equation gives as a solution a great circle.

Surely not an elegant proof like Bar's reference.

But I hope this is correct.

{As a friend of mine pointed out that this set of coordinates is motivated by the fact that the way the "energy" of the curve is being parametrized the z-component of the angular momentum is conserved which is in fact my second Euler-Lagrange equations}