Let $(M, g)$ be a smooth Riemannian manifold and let $\gamma : [0,1] \to M$ be a geodesic. Then in local co-ordinate $(x_1, \dots, x_n)$ where $x_i(t) = \gamma_i(t)$ and we set $v_i(t) = \dot{x}_i(t)$, the goedesic equations are $$ \begin{cases} \dot{x}_k = v_k \\ \dot{v}_k = -\sum_{i,j}^n \Gamma_{ij}^kv_iv_j \end{cases} $$ where $$ \Gamma_{ij}^k = \frac{1}{2}\sum_{\ell = 1}^{n}g^{k\ell}\left(\partial_{x_i}g_{\ell j} + \partial_{x_j}g_{\ell i} - \partial_{x_\ell}g_{ij}\right). $$ are the Christöffel symbols and $(g^{ij})_{ij}$ denotes the inverse matrix of $(g_{ij})_{ij}$.
Now, we identify $TM$ with $T^\ast M$ by setting $\xi_i = \sum_{j=1}^n g_{ij}(x)v_j$. I am trying to show that the geodesic flow on $T^\ast M$ is given by $$ \begin{cases} \dot{x}_k = \partial_{\xi_k}p\\ \dot{\xi}_k = -\partial_{x_k}p. \end{cases} $$ where $p(x,\xi) = \frac12|\xi|_{g(x)}^2$. Assuming the geodesic equations on $TM$, I was able to show that $\dot{x}_k = \partial_{\xi_k}p$.
On the other hand, I am unable to show that $\dot{\xi_k} = -\partial_{x_k}p$. Instead, I obtain $\dot{\xi_k} = \partial_{x_k}p$. Where am I going wrong? My attempt is presented below;
Note that $$ p(x,\xi) = \frac{1}{2}|\xi|^2_{g(x)} = \frac{1}{2}|v|^2_{g(x)} = \frac{1}{2}\sum_{i,j}g_{ij}(x)v_iv_j. $$ Thus, we see that $$ \partial_{x_k}p(x,\xi) = \frac{1}{2}\sum_{i,j}\partial_{x_k}g_{ij}(x)v_iv_j. $$ On the other hand, we have \begin{align*} \dot{\xi}_k &= \sum_{\ell=1}^ng_{k\ell}\dot{v}_\ell + \sum_{i,j}\partial_{x_i}g_{kj}v_i v_j\\ &=-\sum_{i,j,\ell}g_{k\ell}\Gamma^{\ell}_{ij}v_iv_j + \sum_{i,j}\partial_{x_i}g_{kj}v_i v_j\\ &=-\frac{1}{2}\sum_{i,j}\left(\partial_{x_i}g_{kj} + \partial_{x_j}g_{ki} - \partial_{x_k}g_{ij}\right)v_iv_j + \sum_{i,j}\partial_{x_i}g_{kj}v_i v_j\\ &=\frac{1}{2}\sum_{i,j}\left(\partial_{x_i}g_{kj} - \partial_{x_j}g_{ki} +\partial_{x_k}g_{ij}\right)v_iv_j\\ &=\frac{1}{2}\sum_{i,j}\partial_{x_k}g_{ij}v_iv_j = \partial_{x_k}p(x,\xi) \end{align*} for all $k=1,\dots, n$.
Partial derivative differs from full derivative in the way that the result of $\partial/\partial x_1$ is dependent on the choice of basis $x_1,\ldots,x_n$. When you change the basis, the derivative also changes (it forms covector space dual to vector space of coordinates, so it is no surprise). In Lagrange representation, the independent variables are is $x_k$, $v_k$. So when you take $\partial/\partial x_k$, all $v_k$ are assumed constant.
However in Hamilton representation, the independent variables are $x_k$, $\xi_k$. With that being said: $$ 2\partial_{x_k}p=\frac{\partial}{\partial x_k}(g^{ij}\xi_i\xi_j)=\frac{\partial g^{ij}}{\partial x_k}\xi_i\xi_j. $$
The derivative can be found from: $$ 0=\frac{\partial\delta_l^j}{\partial x_k} g^{il}=\frac{\partial(g_{lm}g^{mj})}{\partial x_k} g^{il}= \frac{\partial g_{lm}}{\partial x_k}g^{mj}g^{il} + \frac{\partial g^{mj}}{\partial x_k}g_{lm}g^{il}=\frac{\partial g_{lm}}{\partial x_k}g^{mj}g^{il} + \frac{\partial g^{ij}}{\partial x_k}. $$
So as you can see: $$ 2\left.\frac{\partial p}{\partial x_k}\right|_{\xi_k=\mathop{\mathrm{const}}}=\frac{\partial g^{ij}}{\partial x_k}\xi_i\xi_j = -\frac{\partial g_{lm}}{\partial x_k}g^{mj}g^{il}\xi_i\xi_j= -\frac{\partial g_{lm}}{\partial x_k}v_lv_m $$