I understand that principal curvature https://en.wikipedia.org/wiki/Principal_curvature comes from the eigenvalues of the matrix made of the 2nd derivatives.
Given the two (orthogonal) surfaces tangents $X_u, X_v$ (a surface defined over uv parametric space), for example: $$ X_u= \frac{\partial X}{\partial u}= \begin{bmatrix} x(u) \\ y(u) \\ z(u) \end{bmatrix} \;\;\;\;\;e.g. \begin{bmatrix} u^2-v^2 \\ 2u+2v \\ 3u-v^2 \end{bmatrix} $$
$$ X_v= \frac{\partial X}{\partial v}= \begin{bmatrix} a(v) \\ b(v) \\ c(v) \end{bmatrix}\;\;\;\;\;e.g. \begin{bmatrix} u^2+v^2 \\ v-u^2 \\ v^2-u \end{bmatrix} $$ I can easily find the 2nd derivatives: $X_{uu}, X_{uv},X_{vu},X_{vv}$, e.g. $$ X_{uu}= \begin{bmatrix} 2u\\ 2 \\ 3 \end{bmatrix} $$ but how to write the Hessian matrix for principal curvature and principal direction?
First of all, your examples of $X_u$ and $X_v$ can't occur (note that $X_{uv}\ne X_{vu}$).
You need more information than just the second partial derivatives to get the principal directions and curvatures. You need to calculate the second fundamental form matrix $$\mathbf{II} = \begin{bmatrix} X_{uu}\cdot n & X_{uv}\cdot n \\ X_{vu}\cdot n & X_{vv}\cdot n \end{bmatrix},$$ where $n$ is the unit normal vector, the first fundamental form matrix $$\mathbf{I} = \begin{bmatrix} X_u\cdot X_u & X_u\cdot X_v \\ X_v\cdot X_u & X_v\cdot X_v\end{bmatrix},$$ and then calculate the matrix of the shape operator by taking $\mathbf I^{-1}\mathbf{II}$. This is the matrix whose eigenvalues and eigenvectors you want to find.
By the way, you may find my freely downloadable differential geometry text (see my profile) of use.