Can we symbolically derive geodesic equation for a given metric using its Christoffel symbols?

edit

You can use the method system of class IntegratedGeodesic. For instance, the following code displays the geodesic equation for the round metric of the 2-sphere in polar coordinates $(\theta,\phi)$:

sage: M = manifolds.Sphere(2)
sage: p = M.point((pi/2, pi))
sage: v0 = M.vector(p, (1, 1))
sage: lbda = var('lbda', latex_name=r'\lambda')
sage: C = M.integrated_geodesic(M.metric(), (lbda, 0, 1), v0)
sage: C.system()[0]
[Dphi^2*cos(theta)*sin(theta), -2*Dphi*Dtheta*cos(theta)/sin(theta)]

The two items in the output list are the values of $\ddot\theta$ and $\ddot\phi$, with Dtheta standing for $\dot\theta$ and Dphi standing for $\dot\phi$ and the dot denotes the derivative with respect to the geodesic parameter $\lambda$.