def METHOD_NAME(X1, X2, kernel):
    r"""Computes the matrix of pairwise kernel values for two given datasets.
    Args:
        X1 (list[datapoint]): List of datapoints (first argument)
        X2 (list[datapoint]): List of datapoints (second argument)
        kernel ((datapoint, datapoint) -> float): Kernel function that maps datapoints to kernel value.
    Returns:
        array[float]: The matrix of kernel values.
    **Example:**
    Consider a simple kernel function based on :class:`~.templates.embeddings.AngleEmbedding`:
    .. code-block :: python
        dev = qml.device('default.qubit', wires=2, shots=None)
        @qml.qnode(dev)
        def circuit(x1, x2):
            qml.templates.AngleEmbedding(x1, wires=dev.wires)
            qml.adjoint(qml.templates.AngleEmbedding)(x2, wires=dev.wires)
            return qml.probs(wires=dev.wires)
        kernel = lambda x1, x2: circuit(x1, x2)[0]
    With this method we can systematically evaluate the kernel function ``kernel`` on
    pairs of datapoints, where the points stem from different datasets, like a training
    and a test dataset.
    >>> X_train = np.random.random((4,2))
    >>> X_test = np.random.random((3,2))
    >>> qml.kernels.kernel_matrix(X_train, X_test, kernel)
    tensor([[0.88875298, 0.90655175, 0.89926447],
            [0.93762197, 0.98163781, 0.93076383],
            [0.91977339, 0.9799841 , 0.91582698],
            [0.80376818, 0.98720925, 0.79349212]], requires_grad=True)
    As we can see, for :math:`n` and :math:`m` datapoints in the first and second
    dataset respectively, the output matrix has the shape :math:`n\times m`.
    """
    N = qml.math.shape(X1)[0]
    M = qml.math.shape(X2)[0]
    matrix = qml.math.stack([kernel(x, y) for x, y in product(X1, X2)])
    if qml.math.ndim(matrix[0]) == 0:
        return qml.math.reshape(matrix, (N, M))
    return qml.math.moveaxis(qml.math.reshape(matrix, (N, M, qml.math.size(matrix[0]))), -1, 0)