I am wondering how to actually determine the gradient of a vector in cylindrical coordinates. I have seen a lot of websites that just say what the general form is but I cannot seem to understand how they got there.
The vector in cylindrical coordinates that I am going to use so everyone can follow along is going to be $\vec{V}=V_{r}\hat{r}+V_{\theta}\hat{\theta}+V_{z}\hat{z}$
On any Riemannian manifold (not necessarily curved), the gradient of a function is the metric dual of the exterior derivative.
The exterior derivative relative to any coordinate system of a function is just
$$ \mathrm{d}f = \partial_{x^1} f \mathrm{d}x^1 + \partial_{x^2} f \mathrm{d}x^2 + \cdots + \partial_{x^k} f \mathrm{d}x^k $$
What we need to do, then, is to "hit it with the inverse metric".
In cylindrical coordinates, the metric is
$$ \mathrm{d}r^2 + r^2 \mathrm{d}\theta^2 + \mathrm{d}z^2 $$
which we can write as the matrix $\mathrm{diag}(1, r^2, 1)$. Inverting the matrix gives $\mathrm{diag}(1, r^{-2}, 1)$ and so the inverse metric is
$$ \hat{r}^2 + r^{-2} \hat{\theta}^2 + \hat{z}^2 $$
So applying the inverse metric to the differential form $\mathrm{d}f$ we get
$$ \nabla f = \partial_r f \hat{r} + r^{-2} \partial_\theta f \hat{\theta} + \partial_z f \hat{z} $$
the coefficients I note are computed by taking the inverse metric as a matrix, and multiplying to it the exterior derivative as a column vector:
$$ \begin{pmatrix} 1 \\ & r^{-2} \\ & & 1\end{pmatrix} \begin{pmatrix} \partial_r f \\ \partial_\theta f \\ \partial_z f \end{pmatrix} $$