Somehow I am confused about this.
Say I start with Spherical coordinate $(r,\theta,\phi)$, and I want to find expression $\hat{\phi}$ in terms of $\{\hat{x}, \hat{y}, \hat{z}\}$. At first I thought this is done using Jacobian matrix, but this doesn't seem to be the case. How do I find transformation matrix that turns one basis vectors into another, given the formula $r,\theta,\phi$ in $x,y,z$?
Thank you.
I might be mistaken but I think you are confusing two different uses of the word coordinates. Let me elaborate.
Linear coordinates:
A vector space $V$ may have a basis. For example, a finite dimensional real vector space might have a basis $e_1, \dots, e_n$. Then every vector in $V$ will have coordinates in this basis. Let $v\in V$ and let the coordinates of $v$ be the $n$-tuple $(a_1, \dots, a_n)$ where $a_i \in \mathbb R$. Then $v$ is the vector $a_1 e_1 + \dots + a_n e_n$.
Coordinate charts:
Moving away from linear spaces you can view $\mathbb R^n$ as an $n$-dimensional manifold. An $n$-manifold is a topological space that is locally homeomorphic to open subsets of $\mathbb R^n$. The local homemomorphisms are called charts.
If $n=3$ a chart from the open ball of radius ${1\over 2}$ centered at $(1,1,1)$ in $\mathbb R^3$ is given by the map $f: \mathbb R^3 \to U\subseteq \mathbb R^3$,
$$ (x,y,z) \mapsto (r = \sqrt{x^2 + y^2 + z^2}, \theta = \arccos {z \over \sqrt{x^2 + y^2 + z^2}}, \varphi = \arctan {y \over x})$$
We call $r, \theta$ and $\varphi$ local coordinates.
But the map $f$ above is not linear so it can't be represented by a matrix. Also, it does not represent a change of basis since there is no basis (in the vector space sense) when we talk about charts.