diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py
index 731c43cf..2bfd688f 100644
--- a/pyerrors/__init__.py
+++ b/pyerrors/__init__.py
@@ -487,5 +487,6 @@ from . import linalg
 from . import mpm
 from . import roots
 from . import integrate
+from . import special
 
 from .version import __version__
diff --git a/pyerrors/misc.py b/pyerrors/misc.py
index d71a2bb5..94a7d4c2 100644
--- a/pyerrors/misc.py
+++ b/pyerrors/misc.py
@@ -20,7 +20,7 @@ def print_config():
               "pandas": pd.__version__}
 
     for key, value in config.items():
-        print(f"{key : <10}\t {value}")
+        print(f"{key: <10}\t {value}")
 
 
 def errorbar(x, y, axes=plt, **kwargs):
diff --git a/pyerrors/special.py b/pyerrors/special.py
new file mode 100644
index 00000000..8b36e056
--- /dev/null
+++ b/pyerrors/special.py
@@ -0,0 +1,23 @@
+import scipy
+import numpy as np
+from autograd.extend import primitive, defvjp
+from autograd.scipy.special import j0, y0, j1, y1, jn, yn, i0, i1, iv, ive, beta, betainc, betaln
+from autograd.scipy.special import polygamma, psi, digamma, gamma, gammaln, gammainc, gammaincc, gammasgn, rgamma, multigammaln
+from autograd.scipy.special import erf, erfc, erfinv, erfcinv, logit, expit, logsumexp
+
+
+__all__ = ["beta", "betainc", "betaln",
+           "polygamma", "psi", "digamma", "gamma", "gammaln", "gammainc", "gammaincc", "gammasgn", "rgamma", "multigammaln",
+           "kn", "j0", "y0", "j1", "y1", "jn", "yn", "i0", "i1", "iv", "ive",
+           "erf", "erfc", "erfinv", "erfcinv", "logit", "expit", "logsumexp"]
+
+
+@primitive
+def kn(n, x):
+    """Modified Bessel function of the second kind of integer order n"""
+    if int(n) != n:
+        raise TypeError("The order 'n' needs to be an integer.")
+    return scipy.special.kn(n, x)
+
+
+defvjp(kn, None, lambda ans, n, x: lambda g: - g * 0.5 * (kn(np.abs(n - 1), x) + kn(n + 1, x)))
diff --git a/tests/special_test.py b/tests/special_test.py
new file mode 100644
index 00000000..548e61db
--- /dev/null
+++ b/tests/special_test.py
@@ -0,0 +1,12 @@
+import numpy as np
+import scipy
+import pyerrors as pe
+import pytest
+
+from autograd import jacobian
+from numdifftools import Jacobian as num_jacobian
+
+def test_kn():
+    for n in np.arange(0, 10):
+        for val in np.linspace(0.1, 7.3, 10):
+            assert np.isclose(num_jacobian(lambda x: scipy.special.kn(n, x))(val), jacobian(lambda x: pe.special.kn(n, x))(val), rtol=1e-10, atol=1e-10)