geodesics equation in Riemannian manifold Let $(M,g)$ be a Riemannian metric and $U \subset M$. We know that in local coordinates on $U \subset M$ the equation for a curve to be a geodesic is:

$$
0=γ ̈
=(d^2 γ^k)/(dt^2 ) ∂_k+(dγ^i)/dt   (dγ^j)/dt Γ_{ij}^k ∂_k
$$

Where $Γ_{ij}^k$ are the Christoffel symbols. This is a second-order system of differential equations.

Now are the geodesics always smooth in a Riemannian metric? 
 Thanks in advance.
Note: The Einstein summation convention is in force throughout this answer.End of Note.

Our OP walaa's question, 

Now are the geodesics always smooth in a Riemannian metric?

appears to be somewhat of a misstatement, insofar as smoothness is technically defined with respect to local coordinate patches of a given atlas; both the smoothness of the metric tensor $g_{ij}$ and of the Christoffel symbols $\Gamma_{ij}^k$ are then affirmed if these functions are smooth in such coordinates.  A more carefully stated question along these lines might read:

Are the geodesics of a smooth Riemannian metric smooth?

This is indeed the case, as may be seen by observing that the coefficients $\Gamma_{ij}^k$ occurring in the geodesic equation

$\ddot \gamma = \dfrac{d^2 \gamma^k}{dt^2} \partial_k + \Gamma_{ij}^k \dfrac{d\gamma^i}{dt} \dfrac{d\gamma^j}{dt} \partial_k = 0 \tag 1$

are themselves smooth functions on $M$, being given in terms of the $g_{ij}$ by

$\Gamma_{ij}^k = \dfrac{1}{2}g^{km}(g_{mi, j} + g_{mj, i} - g_{ij, m}), \tag 2$

where

$[g^{ij}] = [g_{ij}]^{-1}; \tag{2.5}$

in these equations we understand that

$\partial_k = \dfrac{\partial}{\partial x_k}, \tag 3$

and

$g_{ij, m} = \dfrac{\partial g_{ij}}{\partial x_m}, \tag 4$

and so forth.

We may cast (1) in the familiar first order form by setting

$\dfrac{d\gamma^i}{dt} = \beta^i, \; 1 \le i \le \dim M; \tag 5$

then

$\dfrac{d^2\gamma^i}{dt^2} = \dfrac{d\beta^i}{dt}; \tag 6$

(1) may then be written

$\dfrac{d\beta^k}{dt}\partial_k + \Gamma_{ij}^k \beta^i \beta^j \partial_k = 0, \tag 7$

which corresponds to the collection of $\dim M$ first order, non-linear ordinary differential equations

$\dfrac{d\beta^k}{dt} + \Gamma_{ij}^k \beta^i \beta^j \ = 0, \; 1 \le k \le \dim M, \tag 8$

that is,

$\dfrac{d\beta^k}{dt} = -\Gamma_{ij}^k \beta^i \beta^j \ = 0, \; 1 \le k \le \dim M; \tag 9$

since the $\Gamma_{ij}^k$ are functions of the position coordinates $x = (x_1, x_2, \ldots, x_{\dim M})$ in $M$, along $\gamma(t)$ we have

$\Gamma_{ij}^k(x) = \Gamma_{ij}^k(\gamma(t)) = \Gamma_{ij}^k(\gamma^1(t), \gamma^2(t), \ldots, \gamma^{\dim M}(t));\tag{10}$

(9) then becomes

$\dfrac{d\beta^k}{dt} = - \Gamma_{ij}^k(\gamma^1(t), \gamma^2(t), \ldots, \gamma^{\dim M}(t))\beta^i \beta^j = 0, \; 1 \le k \le \dim M; \tag{11}$

these together with (5) form a system of $2\dim M$ ordinary, non-linear differential equations for the $\gamma^i(t)$, $\beta^i(t)$, $1 \le i \le \dim M$.

We observe that if the functions $\Gamma_{ij}^k(\gamma^1, \gamma^2, \ldots, \gamma^{\dim M})\beta^i \beta^j$ occurring on the right-hand side of (11) are in fact $C^m$ in the $\gamma^i$, $\beta^i$, then
(5), (11) forms a $C^m$ system of non-linear, ordinary differential equations for the $\gamma^i$, $\beta^i$.  It then follows from the standard theorems that the solution functions $\gamma^i(t)$, $\beta^i(t)$, $1 \le i \le \dim M$ are themselves at least $C^m$;  since the $\Gamma_{ij}^k$ in fact depend on the first derivatives of the $g_{ij}$ (cf. (2)), we may infer that for $C^m$-smooth $g_{ij}$ the $\Gamma_{ij}^k$ are $C^{m - 1}$-smooth, and hence so the $\gamma_i(t)$, $\beta_i(t)$ are at least $C^{m - 1}$.  Since these assertions bind for any $m \ge 1$ it follows that if the metric $g_{ij}$ is $C^\infty$, so are the functions $\gamma_i(t)$, $\beta_i(t)$, $1 \le i \le \dim M$.  Thus the geodesics of a smooth Riemannian metric are themselves smooth.
The geodesic equation says that, locally, geodesics are at least $C^2$. Since the geodesic equation is invariant under coordinates changes, we find out that a geodesic is at least $C^2$ in all of its range. There is no smoothness $C^\infty$ requirement here.