From what I can tell, the error arises because you need to reparametrize your $u_{1}$ curve to have constant speed. Once you do this, the geodesic equations should work out.
Denote the metric by $Edu_{1}du_{1} + Gdu_{2}du_{2}$ (i.e. $E = g_{11}$ and $G = g_{22}$ with your notation above) and observe that $E$ and $G$ depend on only $u_{1}$. As you have expressed above, the geodesic equations reduce to
\begin{align*}
u_{1}^{\prime\prime} + \Gamma^{1}_{11}\left(u^{\prime}\right)^{2} + \Gamma^{1}_{22}\left(v^{\prime}\right)^{2} &= 0\\
u_{2}^{\prime\prime} + \Gamma^{2}_{1 2}u_{1}^{\prime}u_{2}^{\prime} &=0\\
\end{align*}
Since the class of curves under consideration is the meridians, you have correctly observed that all terms involving a factor of $u_{2}^{\prime}$ will vanish. The remaining geodesic equation(s) is
\begin{equation}
u_{1}^{\prime\prime} + \Gamma^{1}_{11} \left(u^{\prime}\right)^{2} = 0,
\end{equation}
where $ \Gamma^{1}_{11} = \frac{f^{\prime}f^{\prime \prime}}{1 + (f^\prime)^2} = \frac{E_{u_{1}}}{{2E}}$
Now given a curve curve of the form $\gamma(t) = \sigma\left(t, c\right)$ (i.e. $u_{1}(t) = t$), we should first reparametrize $\gamma$ to have constant speed. Observe that the arc length function is
$$s(t) = \int\limits_{0}^{t} \vert \vert \gamma^{\prime}(w) \vert \vert dw = \int \limits_{0}^{t} \sqrt{E(w)} dw,$$
where the last equality follows from the fact that the tangent vector of our meridian curve is tangent to the $u_{1}$ curves. Of course, the fundamental theorem of calculus implies that $\frac{ds}{dt} = \sqrt{E(t)} > 0$, and your meridian admits a unit speed reparametrization via the inverse function $t = t(s)$.
Take the unit speed reparametrization to be $\alpha(s) = \sigma(t(s), c)$, (i.e. $u_{1}(s) = u_{1}(t(s))$ where $u_{1}(s) = t(s)$ and note that by the chain rule (and the fact that you have a meridian so you are tangent to a $u_{1}$ curve) you obtain the following derivatives:
$$\frac{du_{1}}{ds} = \frac{du_{1}}{dt}\frac{dt}{ds} = 1\cdot\frac{1}{\sqrt{E}},$$
and
$$\frac{d^{2}u_{1}}{ds^2}= -\frac{E_{u_{1}}}{2E^{2}}.$$
You should now be able to observe that the the unit speed parametrization of your meridian satisfies the geodesic equation above.
Finally, a quick remark regarding the importance of the parametrization. Take the Euclidean plane with metric $ds^{2} = dx^{2} + dy^{2}$ (i.e. with $(g_{ij}) = \begin{pmatrix}
1&0\\
0&1\\
\end{pmatrix}$). If you write down the geodesic equations with the Christoffel symbols, you end up with
$$ x^{\prime \prime} = 0 \hspace{.5in} \textrm{and} \hspace{.5in} y^{\prime \prime} = 0.$$
Now consider the line parametrized by $\gamma(t) = (t^{3}, t^{3})$. The curve $\gamma$ clearly traces out a geodesic, but the given parametrization does not satisfy the geodesic equations.
Key idea. You may simplify your target by choosing an appropriate coordinate before calculations.
Preliminaries. You might know that, for a set of ordinary differential equation (ODE) system, it would not be complete unless you specify initial conditions for them (of course, you can still seek for general solutions even without initial conditions, but this does not help especially if you hope to get a proof with mathematical rigor). So for your ODE system, a complete expression would be
\begin{align}
\ddot{\theta}-\dot{\phi}^2\sin\theta\cos\theta&=0,\\
\ddot{\phi}+2\dot{\theta}\dot{\phi}\cot\theta&=0,\\
\theta|_{\tau=0}&=\theta_0,\\
\dot{\theta}|_{\tau=0}&=u_0,\\
\phi|_{\tau=0}&=\tau_0,\\
\dot{\phi}|_{\tau=0}&=v_0,
\end{align}
where $\theta_0$ and $\phi_0$ are initial values for $\theta$ and $\phi$, while $u_0$ and $v_0$ are initial values for $\dot{\theta}$ and $\dot{\phi}$.
Geometrical interpretation. The initial value $\left(\theta_0,\phi_0\right)$ indicates the point on the sphere at which the mass point starts moving when $\tau=0$, while $\left(u_0,v_0\right)$, a vector pinned to this point, is the velocity of the mass point at $\tau=0$. So $\left(\theta_0,\phi_0,u_0,v_0\right)$ as a whole gives you a tangential vector somewhere on the sphere.
Now recall the spherical coordinate system, as is shown in the following figure. The tricky part is, we can always reconsider the coordinate after the tuple $\left(\theta_0,\phi_0,u_0,v_0\right)$ is given. That is, as $\left(\theta_0,\phi_0,u_0,v_0\right)$ is specified, we perform a transformation, making this $\left(\theta_0,\phi_0,u_0,v_0\right)$ appear to be a vector at the north pole pointing parallel to the $x$-axis in our new coordinate system. Therefore, in our new system, $\theta_0=0$, $\phi_0=0$, $v_0=0$. Plus, to keep the normality of the velocity, let $u_0=1$.
Consequently, our system becomes
\begin{align}
\ddot{\theta}-\dot{\phi}^2\sin\theta\cos\theta&=0,\\
\ddot{\phi}+2\dot{\theta}\dot{\phi}\cot\theta&=0,\\
\theta|_{\tau=0}&=0,\\
\dot{\theta}|_{\tau=0}&=1,\\
\phi|_{\tau=0}&=0,\\
\dot{\phi}|_{\tau=0}&=0.
\end{align}
In addition, to show that the solution corresponds to a great circle, it suffices to show that the above system implies $\phi\equiv 0$ (or more rigorously, $\phi\equiv 0$ when $\theta\in\left[0,\pi\right]$, and $\phi\equiv\pi$ when $\theta\in\left(\pi,2\pi\right)$).
Remark. The symbol "$\equiv$" I adopted here means "always equal to" or "constantly equal to".
Main steps. To begin with, thanks to
\begin{align}
\ddot{\theta}-\dot{\phi}^2\sin\theta\cos\theta&=0,\\
\ddot{\phi}+2\dot{\theta}\dot{\phi}\cot\theta&=0,
\end{align}
we have
$$
\frac{\rm d}{{\rm d}\tau}\left(\dot{\theta}^2+\dot{\phi}^2\sin^2\theta\right)=0,
$$
and the initial values yield
$$
\dot{\theta}^2+\dot{\phi}^2\sin^2\theta\equiv\left(\dot{\theta}^2+\dot{\phi}^2\sin^2\theta\right)\Big|_{\tau=0}=1.
$$
(This is exactly your arc-length-parametrization argument. I mention it here just to highlight that it is a natural consequence of the geodesic equations, for which there is no need to impose it as an additional constraint.)
We observe that, due to $\dot{\theta}|_{\tau=0}=1$, $\theta\not\equiv 0$. Hence $\sin\theta\not\equiv 0$. Thanks to this fact, use the last constraint to eliminate $\dot{\phi}^2$ in the first equation, which leads to
$$
\ddot{\theta}=\left(1-\dot{\theta}^2\right)\cot\theta.
$$
This is a decoupled ODE for $\theta$, and it could be dealt with using the following standard trick. We have
$$
\ddot{\theta}=\frac{\rm d}{{\rm d}\tau}\dot{\theta}=\frac{{\rm d}\theta}{{\rm d}\tau}\frac{{\rm d}\dot{\theta}}{{\rm d}\theta}=\dot{\theta}\frac{{\rm d}\dot{\theta}}{{\rm d}\theta}=\frac{\rm d}{{\rm d}\theta}\left(\frac{1}{2}\dot{\theta}^2\right).
$$
As a result, the decoupled ODE becomes
$$
\frac{\rm d}{{\rm d}\theta}\left(\frac{1}{2}\dot{\theta}^2\right)=\left(1-\dot{\theta}^2\right)\cot\theta\iff-\frac{{\rm d}\left(1-\dot{\theta}^2\right)}{1-\dot{\theta}^2}=2\frac{{\rm d}\sin\theta}{\sin\theta}.
$$
This is nothing but
$$
{\rm d}\log\left[\left(1-\dot{\theta}^2\right)\sin^2\theta\right]=0,
$$
and we immediately have
$$
\left(1-\dot{\theta}^2\right)\sin^2\theta\equiv\left[\left(1-\dot{\theta}^2\right)\sin^2\theta\right]\Big|_{\tau=0}=0,
$$
meaning that, with $\dot{\theta}|_{\tau=0}=1$,
$$
\text{either}\quad\dot{\theta}\equiv 1\quad\text{or}\quad\sin\theta\equiv 0
$$
must be true. As has been mentioned that $\sin\theta\not\equiv 0$, it is a must that $\dot{\theta}\equiv 1$.
Finally, substitute $\theta=\tau$ (obtained from $\dot{\theta}=1$ as well as $\theta|_{\tau=0}=0$) into the second equation, and we have another decoupled system
\begin{align}
\ddot{\phi}+2\dot{\phi}\cot\tau&=0,\\
\phi|_{\tau=0}&=0,\\
\dot{\phi}|_{\tau=0}&=0.
\end{align}
Obviously, $\dot{\phi}\equiv 0$ solves this system. We need to show that this is the only solution this system admits. Suppose for contradiction that there exists some $\dot{\phi}\not\equiv 0$. As such, the equation could be re-expressed as
$$
\frac{{\rm d}\dot{\phi}}{\dot{\phi}}=-2\frac{{\rm d}\sin\tau}{\sin{\tau}}\iff{\rm d}\log\left(\dot{\phi}\sin^2\tau\right)=0,
$$
i.e.,
$$
\dot{\phi}\sin^2\tau=\left(\dot{\phi}\sin^2\tau\right)\Big|_{\tau=0}=0\Longrightarrow\dot{\phi}\equiv 0,
$$
which contradicts our assumption that $\dot{\phi}\not\equiv 0$. As a consequence, $\dot{\phi}\equiv 0$, and thus $\phi\equiv 0$ as per its initial condition, is the only possible solution.
To sum up, we obtain that
\begin{align}
\theta&=\tau,\\
\phi&=0
\end{align}
is the only solution that satisfies the given geodesic equations as well as its specified initial conditions. Obviously, this solution corresponds to a great circle on the sphere.
Therefore, we may conclude that any curve on the unit sphere that satisfies the geodesic equations must be (part of) a great circle.
Best Answer
The Christoffel symbols you dervied are indeed the correct ones for a spherical coordinate system $(r, \theta, \varphi)$.
If you do the same procedure for a system $(r, \varphi, \theta)$ (in the metric tensor, the entries $(22)$ and $(33)$ are now swapped) you will get the Christoffel symbols as stated on Wolfram Mathworld.
It is simply due to the order of $\theta$ and $\varphi$.