Although probably it will be so trivial, I'm actually stuck on last part in this answer https://math.stackexchange.com/a/1166154/712123
I don't know how to prove that $\nabla_{\partial_i} dx^j = -\Gamma_{ki}^j dx^k$ from previous step $\left( \nabla_{\partial_i} dx^j \right)\partial_k = -\Gamma_{ki}^j$, because if I just contract all with $dx^k $ I end up with:
$$\begin{align} dx^k\left( \nabla_{\partial_i} dx^j \right)\partial_k &= -\Gamma_{ki}^jdx^k\\ \left( \nabla_{\partial_i} dx^j \right)\delta^k_k &= -\Gamma_{ki}^jdx^k\\ \left( \nabla_{\partial_i} dx^j \right)&= - \frac{1}{n} \Gamma_{ki}^jdx^k \end{align}$$
where $n$ is the dimension of the manifold. This is because the presence of Dirac delta due to the duality of $dx^k$ and $\partial_k$.
How can I get the desirable answer? Thanks
What you've written makes no sense. You should be considering $$\nabla_{\partial_i}\big(dx^j\partial_k\big) = \big(\nabla_{\partial_i}dx^j\big)\partial_k + dx^j\big(\nabla_{\partial_i}\partial_k\big).$$ This is then $$\nabla_{\partial_i}\delta^j_k = 0 = \big(\nabla_{\partial_i}dx^j\big)\partial_k + dx^j\big(\Gamma_{ki}^\ell\partial_\ell\big) = \big(\nabla_{\partial_i}dx^j\big)\partial_k + \Gamma_{ki}^j.$$ Therefore, $\nabla_{\partial_i}dx^j = -\Gamma_{ki}^j dx^k$.
EDIT: To get the last step, write $\nabla_{\partial_i}dx^j = a_{i\ell}^j dx^\ell$ for some $a_{i\ell}^j$. Evaluate both sides on $\partial_k$. You get $$-\Gamma_{ki}^j=\big(\nabla_{\partial_i}dx^j\big)\partial_k = a_{i\ell}^j dx^\ell\partial_k = a_{i\ell}^j\delta^\ell_k = a_{ik}^j.$$
EDIT: The mistake in your computation occurs at the first step. There is no "contraction" to do. The left-hand side is a scalar, not a vector. Note that $\nabla_{\partial_i}dx^j = \omega_i^j$ is a $1$-form and you're contracting it with $\partial_k$ to start with. You're trying to say that if we take another $1$-form $\phi$ ($dx^k$ in your case), then $(\omega_i^j\partial_k)\phi = (\phi\partial_k)\omega_i^j$. This is just not remotely valid. Let's just try a simple example: $$(dx^i\partial_k)dx^m \overset{?}= (dx^m\partial_k)dx^i \iff \delta^i_k dx^m \overset{?}= \delta^m_k dx^i.$$ Well, this holds only if $i=m$, or else you have linearly independent $1$-forms equal to one another. To emphasize, you think you're contracting a $1$-form with a vector, but you're not able to do so — you're multiplying a $1$-form by a scalar function, and there's no way to reassociate and contract $dx^k$ with $\partial_k$.