Finite differences
Let's take a general $2n+1$ point central finite difference rule with uniform grid spacing
$$ f'(x) \approx h^{-1} \sum_{i=-n}^{n} w_i f(x+ ih), $$
such as two-point central finite differences with $n=1$ and ${\bf w}=(-1/2, 0, 1/2)$ which yields the approximation
$$ f'(x) \approx \frac {f(x+h)-f(x-h)} {2h}. $$
Geometry optimization
For geometrical gradients, our function is the energy $E(\{\bf R\})$ where ${\bf R}_I = (X_I, Y_I, Z_I)$ are the coordinates of the $I$-th nucleus. (If we have $N$ nuclei, $\{\bf R\}$ is a vector with $3N$ elements, consisting of the $x$, $y$, and $z$ coordinates of each $N$ nuclei.)
The $k$-th component ($x$, $y$ or $z$) of the force acting on nucleus $I$ is
$$ F_{Ik} = -\frac {\partial E} {\partial R_{Ik}}, $$ and it can be evaluated efficiently with analytical methods for a lot of quantum chemical methods. Note that here I use two indices: one for the nucleus and another for the cartesian component.
Simplification of the notation with compound indices
I will now simplify the notation $Ik \to i$. If you number the nuclei from 0 to N-1, then the x, y, and z coordinates of the Ith nucleus will be stored at indices $i=3I$, $i=3I+1$ and $i=3I+2$, respectively. I will also denote the vector of nuclear coordinates with ${\bf x}$ and then the force is
$$ F_{i} = -\frac {\partial E} {\partial x_i}. $$
Hessian
The Hessian is the second derivative, and like you said, it is a $3N \times 3N$ square matrix. Analytical Hessians are less often available than analytical gradients, and they can be evaluated either seminumerically, if analytical gradients are available, or fully numerically by applying two rounds of finite differences.
The Hessian is
$$ H_{ij} = \frac {\partial^2 E} {\partial x_i \partial x_j} $$
where I now use the compound indexing.
Seminumerical scheme
Note that if I know the gradients $g_i = \partial E / \partial x_i = -F_i$, I can evaluate the Hessian with finite differences as
$$ H_{ij} = \frac {\partial^2 E} {\partial x_i \partial x_j} =\frac {\partial g_i} {\partial x_j} $$
This form is explicit, since it tells you that the Hessian measures how much the $i$th gradient changes when the $j$th coordinate is modified. How is this evaluated? Let's take the two-point equation above for simplicity
$$ H_{ij} =\frac {\partial g_i} {\partial x_j} \approx \frac {g_i({\bf x}+h\hat{e}_j)-g_i({\bf x}-h\hat{e}_j)} {2h}. $$
That is, for each column $j$ in the Hessian, you get the row $i$ by taking the difference of gradients of all the atoms where the gradients are evaluated with coordinates where the $j$-th coordinate has been modified. This means that for this rule, you need to do 2 gradient evaluations for every $3N$ columns, or $6N$ gradient calculations.
Note in the above that if the $i$ derivative is evaluated analytically and the $j$ derivative with finite differences, the Hessian matrix will have numerical noise in column indices but not in row indices, and thereby will be asymmetric, even though it is supposed to be symmetric. This is why it is a good idea to re-symmetrize the Hessian, $$ H_{ij} \to \frac {H_{ij} + H_{ji}} 2 $$.
Fully numerical scheme
What about the fully numerical Hessian? It is nothing more complicated than the above; you just use the finite difference equation again to expand $g_i$ to get
$$ H_{ij} \approx \frac {E({\bf x}+h\hat{e}_i+h\hat{e}_j)-E({\bf x}-h\hat{e}_i+h\hat{e}_j)-E({\bf x}+h\hat{e}_i-h\hat{e}_j)+E({\bf x}-h\hat{e}_i-h\hat{e}_j)} {4h^2}. $$
Even though this example only used two-point rules, you can use higher-order ones as well, such as the ones presented in Wikipedia.
