diff --git a/Cargo.lock b/Cargo.lock
index df09935..99604b4 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -43,6 +43,15 @@ dependencies = [
  "libc",
 ]
 
+[[package]]
+name = "approx"
+version = "0.5.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6"
+dependencies = [
+ "num-traits",
+]
+
 [[package]]
 name = "arrayvec"
 version = "0.7.6"
@@ -443,6 +452,7 @@ dependencies = [
  "ryu",
  "serde",
  "serde_json",
+ "statrs",
 ]
 
 [[package]]
@@ -965,6 +975,16 @@ version = "0.3.11"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "38b58827f4464d87d377d175e90bf58eb00fd8716ff0a62f80356b5e61555d0d"
 
+[[package]]
+name = "statrs"
+version = "0.18.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2a3fe7c28c6512e766b0874335db33c94ad7b8f9054228ae1c2abd47ce7d335e"
+dependencies = [
+ "approx",
+ "num-traits",
+]
+
 [[package]]
 name = "subtle"
 version = "2.5.0"
diff --git a/base/Cargo.toml b/base/Cargo.toml
index 7befc11..66433c4 100644
--- a/base/Cargo.toml
+++ b/base/Cargo.toml
@@ -19,6 +19,7 @@ regex = { version = "1.0", optional = true}
 regex-lite = { version = "0.1.6", optional = true}
 bitcode = "0.6.3"
 csv = "1.3.0"
+statrs = { version = "0.18.0", default-features = false, features = [] }
 
 [features]
 default = ["use_regex_full"]
diff --git a/base/src/functions/engineering/bessel.rs b/base/src/functions/engineering/bessel.rs
index a461abd..12a09a1 100644
--- a/base/src/functions/engineering/bessel.rs
+++ b/base/src/functions/engineering/bessel.rs
@@ -1,10 +1,12 @@
+use statrs::function::erf::erf;
+
 use crate::{
     calc_result::CalcResult,
     expressions::{parser::Node, token::Error, types::CellReferenceIndex},
     model::Model,
 };
 
-use super::transcendental::{bessel_i, bessel_j, bessel_k, bessel_y, erf};
+use super::transcendental::{bessel_i, bessel_j, bessel_k, bessel_y};
 // https://root.cern/doc/v610/TMath_8cxx_source.html
 
 // Notice that the parameters for Bessel functions in Excel and here have inverted order
diff --git a/base/src/functions/engineering/transcendental/erf.rs b/base/src/functions/engineering/transcendental/erf.rs
deleted file mode 100644
index 9c56a4c..0000000
--- a/base/src/functions/engineering/transcendental/erf.rs
+++ /dev/null
@@ -1,53 +0,0 @@
-pub(crate) fn erf(x: f64) -> f64 {
-    let cof = vec![
-        -1.3026537197817094,
-        6.419_697_923_564_902e-1,
-        1.9476473204185836e-2,
-        -9.561_514_786_808_63e-3,
-        -9.46595344482036e-4,
-        3.66839497852761e-4,
-        4.2523324806907e-5,
-        -2.0278578112534e-5,
-        -1.624290004647e-6,
-        1.303655835580e-6,
-        1.5626441722e-8,
-        -8.5238095915e-8,
-        6.529054439e-9,
-        5.059343495e-9,
-        -9.91364156e-10,
-        -2.27365122e-10,
-        9.6467911e-11,
-        2.394038e-12,
-        -6.886027e-12,
-        8.94487e-13,
-        3.13092e-13,
-        -1.12708e-13,
-        3.81e-16,
-        7.106e-15,
-        -1.523e-15,
-        -9.4e-17,
-        1.21e-16,
-        -2.8e-17,
-    ];
-
-    let mut d = 0.0;
-    let mut dd = 0.0;
-
-    let x_abs = x.abs();
-
-    let t = 2.0 / (2.0 + x_abs);
-    let ty = 4.0 * t - 2.0;
-
-    for j in (1..=cof.len() - 1).rev() {
-        let tmp = d;
-        d = ty * d - dd + cof[j];
-        dd = tmp;
-    }
-
-    let res = t * f64::exp(-x_abs * x_abs + 0.5 * (cof[0] + ty * d) - dd);
-    if x < 0.0 {
-        res - 1.0
-    } else {
-        1.0 - res
-    }
-}
diff --git a/base/src/functions/engineering/transcendental/mod.rs b/base/src/functions/engineering/transcendental/mod.rs
index 9a90a48..ea01fbe 100644
--- a/base/src/functions/engineering/transcendental/mod.rs
+++ b/base/src/functions/engineering/transcendental/mod.rs
@@ -4,7 +4,6 @@ mod bessel_j1_y1;
 mod bessel_jn_yn;
 mod bessel_k;
 mod bessel_util;
-mod erf;
 
 #[cfg(test)]
 mod test_bessel;
@@ -13,4 +12,3 @@ pub(crate) use bessel_i::bessel_i;
 pub(crate) use bessel_jn_yn::jn as bessel_j;
 pub(crate) use bessel_jn_yn::yn as bessel_y;
 pub(crate) use bessel_k::bessel_k;
-pub(crate) use erf::erf;