code: def kernel_matrix(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,
      )
  
label_desc: Use image node_id 2 for calling a global function with example usage: kernel_matrix(X1, X2, kernel) and returns: qml, qml