diff --git a/CHANGELOG.md b/CHANGELOG.md
index f2ad67edf49..1860bb8aafa 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -19,6 +19,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * Added implementation of `dpnp.piecewise` [#2550](https://github.com/IntelPython/dpnp/pull/2550)
 * Added implementation of `dpnp.linalg.lu_solve` for 2D inputs (SciPy-compatible) [#2575](https://github.com/IntelPython/dpnp/pull/2575)
 * Added implementation of `dpnp.special.erfc` [#2588](https://github.com/IntelPython/dpnp/pull/2588)
+* Added implementation of `dpnp.special.erfcx` [#2596](https://github.com/IntelPython/dpnp/pull/2596)
 
 ### Changed
 
diff --git a/doc/reference/special.rst b/doc/reference/special.rst
index 5e6735a6afb..3efc836fa89 100644
--- a/doc/reference/special.rst
+++ b/doc/reference/special.rst
@@ -15,6 +15,5 @@ Error function and Fresnel integrals
    erf
    erfc
    erfcx
-   erfi
    erfinv
    erfcinv
diff --git a/dpnp/backend/extensions/ufunc/elementwise_functions/erf_funcs.cpp b/dpnp/backend/extensions/ufunc/elementwise_functions/erf_funcs.cpp
index 1a441d23e83..1f6f15be424 100644
--- a/dpnp/backend/extensions/ufunc/elementwise_functions/erf_funcs.cpp
+++ b/dpnp/backend/extensions/ufunc/elementwise_functions/erf_funcs.cpp
@@ -209,6 +209,7 @@ static void populate(py::module_ m,
 
 MACRO_DEFINE_IMPL(erf, Erf);
 MACRO_DEFINE_IMPL(erfc, Erfc);
+MACRO_DEFINE_IMPL(erfcx, Erfcx);
 } // namespace impl
 
 void init_erf_funcs(py::module_ m)
@@ -228,5 +229,9 @@ void init_erf_funcs(py::module_ m)
     impl::populate<impl::ErfcContigFactory, impl::ErfcStridedFactory>(
         m, "_erfc", "", impl::erfc_contig_dispatch_vector,
         impl::erfc_strided_dispatch_vector);
+
+    impl::populate<impl::ErfcxContigFactory, impl::ErfcxStridedFactory>(
+        m, "_erfcx", "", impl::erfcx_contig_dispatch_vector,
+        impl::erfcx_strided_dispatch_vector);
 }
 } // namespace dpnp::extensions::ufunc
diff --git a/dpnp/backend/extensions/vm/erf_funcs.cpp b/dpnp/backend/extensions/vm/erf_funcs.cpp
index 2414c56a9e5..4630ccc10be 100644
--- a/dpnp/backend/extensions/vm/erf_funcs.cpp
+++ b/dpnp/backend/extensions/vm/erf_funcs.cpp
@@ -130,6 +130,7 @@ using ew_cmn_ns::unary_contig_impl_fn_ptr_t;
 
 MACRO_DEFINE_IMPL(erf, Erf);
 MACRO_DEFINE_IMPL(erfc, Erfc);
+MACRO_DEFINE_IMPL(erfcx, Erfcx);
 
 template <template <typename fnT, typename T> typename factoryT>
 static void populate(py::module_ m,
@@ -184,5 +185,11 @@ void init_erf_funcs(py::module_ m)
         "Call `erfc` function from OneMKL VM library to compute the "
         "complementary error function value of vector elements",
         impl::erfc_contig_dispatch_vector);
+
+    impl::populate<impl::ErfcxContigFactory>(
+        m, "_erfcx",
+        "Call `erfcx` function from OneMKL VM library to compute the scaled "
+        "complementary error function value of vector elements",
+        impl::erfcx_contig_dispatch_vector);
 }
 } // namespace dpnp::extensions::vm
diff --git a/dpnp/backend/kernels/elementwise_functions/erf.hpp b/dpnp/backend/kernels/elementwise_functions/erf.hpp
index a292a0a3f3f..2b6203580a0 100644
--- a/dpnp/backend/kernels/elementwise_functions/erf.hpp
+++ b/dpnp/backend/kernels/elementwise_functions/erf.hpp
@@ -29,6 +29,19 @@
 
 #include <sycl/sycl.hpp>
 
+/**
+ * Include <sycl/ext/intel/math.hpp> only when targeting to Intel devices.
+ */
+#if defined(__INTEL_LLVM_COMPILER)
+#define __SYCL_EXT_INTEL_MATH_SUPPORT
+#endif
+
+#if defined(__SYCL_EXT_INTEL_MATH_SUPPORT)
+#include <sycl/ext/intel/math.hpp>
+#else
+#include <cmath>
+#endif
+
 namespace dpnp::kernels::erfs
 {
 template <typename OpT, typename ArgT, typename ResT>
@@ -62,13 +75,20 @@ struct BaseFunctor
         template <typename Tp>                                                 \
         static Tp apply(const Tp &x)                                           \
         {                                                                      \
-            return sycl::__name__(x);                                          \
+            return __name__(x);                                                \
         }                                                                      \
     };                                                                         \
                                                                                \
     template <typename ArgT, typename ResT>                                    \
     using __f_name__##Functor = BaseFunctor<__f_name__##Op, ArgT, ResT>;
 
-MACRO_DEFINE_FUNCTOR(erf, Erf);
-MACRO_DEFINE_FUNCTOR(erfc, Erfc);
+MACRO_DEFINE_FUNCTOR(sycl::erf, Erf);
+MACRO_DEFINE_FUNCTOR(sycl::erfc, Erfc);
+MACRO_DEFINE_FUNCTOR(
+#if defined(__SYCL_EXT_INTEL_MATH_SUPPORT)
+    sycl::ext::intel::math::erfcx,
+#else
+    std::erfc,
+#endif
+    Erfcx);
 } // namespace dpnp::kernels::erfs
diff --git a/dpnp/special/__init__.py b/dpnp/special/__init__.py
index 8e1a3fcf24e..a50582e8009 100644
--- a/dpnp/special/__init__.py
+++ b/dpnp/special/__init__.py
@@ -42,9 +42,11 @@
 from ._erf import (
     erf,
     erfc,
+    erfcx,
 )
 
 __all__ = [
     "erf",
     "erfc",
+    "erfcx",
 ]
diff --git a/dpnp/special/_erf.py b/dpnp/special/_erf.py
index 6060f09e943..b7c3623e35e 100644
--- a/dpnp/special/_erf.py
+++ b/dpnp/special/_erf.py
@@ -42,7 +42,7 @@
 import dpnp.backend.extensions.ufunc._ufunc_impl as ufi
 from dpnp.dpnp_algo.dpnp_elementwise_common import DPNPUnaryFunc
 
-__all__ = ["erf", "erfc"]
+__all__ = ["erf", "erfc", "erfcx"]
 
 
 # pylint: disable=too-few-public-methods
@@ -96,7 +96,6 @@ def __call__(self, x, out=None):  # pylint: disable=signature-differs
 :obj:`dpnp.special.erfinv` : Inverse of the error function.
 :obj:`dpnp.special.erfcinv` : Inverse of the complementary error function.
 :obj:`dpnp.special.erfcx` : Scaled complementary error function.
-:obj:`dpnp.special.erfi` : Imaginary error function.
 
 Notes
 -----
@@ -152,7 +151,6 @@ def __call__(self, x, out=None):  # pylint: disable=signature-differs
 :obj:`dpnp.special.erfinv` : Inverse of the error function.
 :obj:`dpnp.special.erfcinv` : Inverse of the complementary error function.
 :obj:`dpnp.special.erfcx` : Scaled complementary error function.
-:obj:`dpnp.special.erfi` : Imaginary error function.
 
 Examples
 --------
@@ -171,3 +169,48 @@ def __call__(self, x, out=None):  # pylint: disable=signature-differs
     mkl_fn_to_call="_mkl_erf_to_call",
     mkl_impl_fn="_erfc",
 )
+
+_ERFCX_DOCSTRING = r"""
+Calculates the scaled complementary error function of a given input array.
+
+It is defined as :math:`\exp(x^2) * \operatorname{erfc}(x)`.
+
+For full documentation refer to :obj:`scipy.special.erfcx`.
+
+Parameters
+----------
+x : {dpnp.ndarray, usm_ndarray}
+    Input array, expected to have a real-valued floating-point data type.
+out : {dpnp.ndarray, usm_ndarray}, optional
+    Optional output array for the function values.
+
+Returns
+-------
+out : dpnp.ndarray
+    The values of the scaled complementary error function at the given points
+    `x`.
+
+See Also
+--------
+:obj:`dpnp.special.erf` : Gauss error function.
+:obj:`dpnp.special.erfc` : Complementary error function.
+:obj:`dpnp.special.erfinv` : Inverse of the error function.
+:obj:`dpnp.special.erfcinv` : Inverse of the complementary error function.
+
+Examples
+--------
+>>> import dpnp as np
+>>> x = np.linspace(-3, 3, num=4)
+>>> np.special.erfcx(x)
+array([1.62059889e+04, 5.00898008e+00, 4.27583576e-01, 1.79001151e-01])
+
+"""
+
+erfcx = DPNPErf(
+    "erfcx",
+    ufi._erf_result_type,
+    ufi._erfcx,
+    _ERFCX_DOCSTRING,
+    mkl_fn_to_call="_mkl_erf_to_call",
+    mkl_impl_fn="_erfcx",
+)
diff --git a/dpnp/tests/test_special.py b/dpnp/tests/test_special.py
index b7eb159b127..462bbd42cd7 100644
--- a/dpnp/tests/test_special.py
+++ b/dpnp/tests/test_special.py
@@ -14,7 +14,7 @@
 
 
 @with_requires("scipy")
-@pytest.mark.parametrize("func", ["erf", "erfc"])
+@pytest.mark.parametrize("func", ["erf", "erfc", "erfcx"])
 class TestCommon:
     @pytest.mark.parametrize(
         "dt", get_all_dtypes(no_none=True, no_float16=False, no_complex=True)
@@ -65,18 +65,36 @@ def test_complex(self, func, dt):
 
 class TestConsistency:
 
-    def test_erfc(self):
+    tol = 8 * dpnp.finfo(dpnp.default_float_type()).resolution
+
+    def _check_variant_func(self, func, other_func, rtol, atol=0):
         # TODO: replace with dpnp.random.RandomState, once pareto is added
         rng = numpy.random.RandomState(1234)
         n = 10000
         a = rng.pareto(0.02, n) * (2 * rng.randint(0, 2, n) - 1)
         a = dpnp.array(a)
+        a = a[::-1]
 
-        res = 1 - dpnp.special.erf(a)
+        res = other_func(a)
         mask = dpnp.isfinite(res)
         a = a[mask]
 
-        tol = 8 * dpnp.finfo(a).resolution
-        assert dpnp.allclose(
-            dpnp.special.erfc(a), res[mask], rtol=tol, atol=tol
+        x, y = func(a), res[mask]
+        if not dpnp.allclose(x, y, rtol=rtol, atol=atol):
+            # calling numpy testing func, because it's more verbose
+            assert_allclose(x.asnumpy(), y.asnumpy(), rtol=rtol, atol=atol)
+
+    def test_erfc(self):
+        self._check_variant_func(
+            dpnp.special.erfc,
+            lambda z: 1 - dpnp.special.erf(z),
+            rtol=self.tol,
+            atol=self.tol,
+        )
+
+    def test_erfcx(self):
+        self._check_variant_func(
+            dpnp.special.erfcx,
+            lambda z: dpnp.exp(z * z) * dpnp.special.erfc(z),
+            rtol=10 * self.tol,
         )
diff --git a/dpnp/tests/test_strides.py b/dpnp/tests/test_strides.py
index 607934c472e..fc874c674dd 100644
--- a/dpnp/tests/test_strides.py
+++ b/dpnp/tests/test_strides.py
@@ -167,7 +167,7 @@ def test_reduce_hypot(dtype, stride):
 
 
 @with_requires("scipy")
-@pytest.mark.parametrize("func", ["erf", "erfc"])
+@pytest.mark.parametrize("func", ["erf", "erfc", "erfcx"])
 @pytest.mark.parametrize("stride", [2, -1, -3])
 def test_erf_funcs(func, stride):
     import scipy.special
diff --git a/dpnp/tests/test_sycl_queue.py b/dpnp/tests/test_sycl_queue.py
index 1ece65860d5..7123a2b9a72 100644
--- a/dpnp/tests/test_sycl_queue.py
+++ b/dpnp/tests/test_sycl_queue.py
@@ -1489,7 +1489,7 @@ def test_interp(device, left, right, period):
     assert_sycl_queue_equal(result.sycl_queue, x.sycl_queue)
 
 
-@pytest.mark.parametrize("func", ["erf", "erfc"])
+@pytest.mark.parametrize("func", ["erf", "erfc", "erfcx"])
 @pytest.mark.parametrize("device", valid_dev, ids=dev_ids)
 def test_erf_funcs(func, device):
     x = dpnp.linspace(-3, 3, num=5, device=device)
diff --git a/dpnp/tests/test_usm_type.py b/dpnp/tests/test_usm_type.py
index e97945e111a..cd644af794a 100644
--- a/dpnp/tests/test_usm_type.py
+++ b/dpnp/tests/test_usm_type.py
@@ -1296,7 +1296,7 @@ def test_choose(usm_type_x, usm_type_ind):
     assert z.usm_type == du.get_coerced_usm_type([usm_type_x, usm_type_ind])
 
 
-@pytest.mark.parametrize("func", ["erf", "erfc"])
+@pytest.mark.parametrize("func", ["erf", "erfc", "erfcx"])
 @pytest.mark.parametrize("usm_type", list_of_usm_types)
 def test_erf_funcs(func, usm_type):
     x = dpnp.linspace(-3, 3, num=5, usm_type=usm_type)
diff --git a/dpnp/tests/third_party/cupyx/scipy_tests/special_tests/test_erf.py b/dpnp/tests/third_party/cupyx/scipy_tests/special_tests/test_erf.py
index c4a640a97d9..5fcc3d3e452 100644
--- a/dpnp/tests/third_party/cupyx/scipy_tests/special_tests/test_erf.py
+++ b/dpnp/tests/third_party/cupyx/scipy_tests/special_tests/test_erf.py
@@ -27,7 +27,6 @@ def test_erf(self):
     def test_erfc(self):
         self.check_unary("erfc")
 
-    @pytest.mark.skip("erfcx() is not supported yet")
     @testing.with_requires("scipy>=1.16.0")
     def test_erfcx(self):
         self.check_unary("erfcx")