diff --git a/doc/whats_new/v1.2.rst b/doc/whats_new/v1.2.rst
index 93081fbfd44b4..4f965ea46be65 100644
--- a/doc/whats_new/v1.2.rst
+++ b/doc/whats_new/v1.2.rst
@@ -198,6 +198,13 @@ Changelog
 :mod:`sklearn.feature_selection`
 ................................
 
+:mod:`sklearn.gaussian_process`
+...............................
+
+- |Fix| Fix :class:`gaussian_process.kernels.Matern` gradient computation with
+  `nu=0.5` for PyPy (and possibly other non CPython interpreters). :pr:`24245`
+  by :user:`Loïc Estève <lesteve>`.
+
 :mod:`sklearn.linear_model`
 ...........................
 
diff --git a/sklearn/gaussian_process/kernels.py b/sklearn/gaussian_process/kernels.py
index 3ac3866cf9ba7..1e0866afb6a4d 100644
--- a/sklearn/gaussian_process/kernels.py
+++ b/sklearn/gaussian_process/kernels.py
@@ -1733,9 +1733,14 @@ def __call__(self, X, Y=None, eval_gradient=False):
 
             if self.nu == 0.5:
                 denominator = np.sqrt(D.sum(axis=2))[:, :, np.newaxis]
-                K_gradient = K[..., np.newaxis] * np.divide(
-                    D, denominator, where=denominator != 0
+                divide_result = np.zeros_like(D)
+                np.divide(
+                    D,
+                    denominator,
+                    out=divide_result,
+                    where=denominator != 0,
                 )
+                K_gradient = K[..., np.newaxis] * divide_result
             elif self.nu == 1.5:
                 K_gradient = 3 * D * np.exp(-np.sqrt(3 * D.sum(-1)))[..., np.newaxis]
             elif self.nu == 2.5:
diff --git a/sklearn/gaussian_process/tests/test_kernels.py b/sklearn/gaussian_process/tests/test_kernels.py
index 05dc682b13e9c..56ab9c8b6c2bf 100644
--- a/sklearn/gaussian_process/tests/test_kernels.py
+++ b/sklearn/gaussian_process/tests/test_kernels.py
@@ -34,7 +34,6 @@
     assert_array_equal,
     assert_array_almost_equal,
     assert_allclose,
-    fails_if_pypy,
 )
 
 
@@ -70,8 +69,6 @@
     kernels.append(PairwiseKernel(gamma=1.0, metric=metric))
 
 
-# Numerical precisions errors in PyPy
-@fails_if_pypy
 @pytest.mark.parametrize("kernel", kernels)
 def test_kernel_gradient(kernel):
     # Compare analytic and numeric gradient of kernels.