It appears that the matrix from rbf is a correlation matrix, with ones on the diagonal. It's also a real-square symmetric matrix, since the upper triangular (off-diagonal elements) are a mirror image of the lower triangular.
The cov matrix looks like a covariance matrix, with diagonal elements that are the variance of each column, and off-diagonals are covariances. You can check if the following are true:
$\sigma_{1}^2 = \frac{\sum_{i=1}^3 (x_{i1} - \bar{x}_1 )^2}{n-1} = 0.1803 \rightarrow\textrm{variance of feature 1}$
$\sigma_{12} = \frac{\sum_{i=1}^3 (x_{i1} - \bar{x}_1 )(x_{i2} - \bar{x}_2 )}{n-1} = 0.4715 \rightarrow\textrm{covariance between feature 1 and 2}$
If so, the cov is the covariance matrix.
I would not commingle kernel methods with the covariance matrix, mostly because kernel tricks could throw you off. If you stick with statistical notation and calculations for obtaining the covariance matrix:
\begin{equation} \boldsymbol{\sigma} = \begin{bmatrix} \sigma_1^2 & \sigma_{12} & \cdots & \sigma_{14} \\ \sigma_{21} & \sigma_2^2 & \cdots & \sigma_{24} \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_{41} & \sigma_{42} & \cdots & \sigma_4^2 \\ \end{bmatrix}, \end{equation}
then you can't go wrong. The correlation matrix is simply:
\begin{equation} \boldsymbol{\rho} = \begin{bmatrix} 1 & \rho_{12} & \cdots & \rho_{14} \\ \rho_{21} & 1 & \cdots & \rho_{24} \\ \vdots & \vdots & \ddots & \vdots \\ \rho_{41} & \rho_{42} & \cdots & 1 \\ \end{bmatrix}, \end{equation}
where all diagonal elements are ones, and elements in the off-diagonal are calculated as
$\rho_{jk} = \frac{\sigma_{jk}}{\sigma_j \sigma_j}$,
where $\sigma_j $ = $\sqrt{\sigma_j^2}$, which are on the diagonal of the covariance matrix.
Now, more to the answer for what you did. You simulated standard normal variates, $z_{ij}$, which have mean zero and variance unity, i.e. $N(0,1)$. The correlation between feature $j$ and feature $k$ (in columns) is defined as
$\rho_{jk} = \frac{\sum_{i=1}^n \left( \frac{x_{ij} - \bar{x}_j}{\sigma_j} \right) \left( \frac{x_{ik} - \bar{x}_k}{\sigma_k} \right) }{n-1} $,
which is simply equal to the product of two z-variates, divided by $n-1$:
$\rho_{jk} = \frac{\sum_i^n z_{ij} z_{ik}}{n-1} $,
since
$z_{ij} = \frac{x_{ij} - \bar{x}_j}{\sigma_j}$.
So your notation should be:
$\boldsymbol{\rho} = \frac{\mathbf{Z}^\top \mathbf{Z}}{n-1}$.
With $n$ rows and $p$ columns in a $\mathbf{Z}$ matrix, plugging in the dimensions, gives
$\underset{p \times p}{\boldsymbol{\rho}} = \frac{\underset{p \times n}{\mathbf{Z}^\top} \underset{n \times p}{\mathbf{Z}}}{n-1}$.
So for column-based correlation matrices, the tranposed matrix, $\mathbf{Z}^\top$, is on the left, and for column-based covariance matrices, $\mathbf{X}^\top$, is also on the left.
It does look like you simulated a $(3 \times 5)$ matrix, and since your matrices are ($3 \times 3)$, you ran row-based correlation and covariance -- which yields ($n \times n)$ matrices. In statistics, it's more common to run correlation and covariance on columns to yield ($p \times p)$ matrices, where $j=1,2,\ldots,p$ represents the features (in columns).