diff --git a/arb_hypgeom.h b/arb_hypgeom.h
index ab992d50..c81b8fec 100644
--- a/arb_hypgeom.h
+++ b/arb_hypgeom.h
@@ -68,6 +68,12 @@ void arb_hypgeom_li(arb_t res, const arb_t z, int offset, slong prec);
 void _arb_hypgeom_li_series(arb_ptr g, arb_srcptr h, slong hlen, int offset, slong len, slong prec);
 void arb_hypgeom_li_series(arb_poly_t g, const arb_poly_t h, int offset, slong len, slong prec);
 
+void arb_hypgeom_bessel_j(arb_t res, const arb_t nu, const arb_t z, slong prec);
+void arb_hypgeom_bessel_y(arb_t res, const arb_t nu, const arb_t z, slong prec);
+void arb_hypgeom_bessel_jy(arb_t res1, arb_t res2, const arb_t nu, const arb_t z, slong prec);
+void arb_hypgeom_bessel_i(arb_t res, const arb_t nu, const arb_t z, slong prec);
+void arb_hypgeom_bessel_k(arb_t res, const arb_t nu, const arb_t z, slong prec);
+
 void arb_hypgeom_airy(arb_t ai, arb_t aip, arb_t bi, arb_t bip, const arb_t z, slong prec);
 void arb_hypgeom_airy_jet(arb_ptr ai, arb_ptr bi, const arb_t z, slong len, slong prec);
 void arb_hypgeom_airy_series(arb_poly_t ai, arb_poly_t ai_prime,
diff --git a/arb_hypgeom/wrappers.c b/arb_hypgeom/wrappers.c
index a22702f6..53dbc5c4 100644
--- a/arb_hypgeom/wrappers.c
+++ b/arb_hypgeom/wrappers.c
@@ -301,3 +301,91 @@ arb_hypgeom_pfq(arb_t res, arb_srcptr a, slong p, arb_srcptr b, slong q, const a
     _acb_vec_clear(t, p + q + 1);
 }
 
+void
+arb_hypgeom_bessel_j(arb_t res, const arb_t nu, const arb_t z, slong prec)
+{
+    acb_t t, u;
+    acb_init(t);
+    acb_init(u);
+    arb_set(acb_realref(t), nu);
+    arb_set(acb_realref(u), z);
+    acb_hypgeom_bessel_j(t, t, u, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+}
+
+void
+arb_hypgeom_bessel_y(arb_t res, const arb_t nu, const arb_t z, slong prec)
+{
+    acb_t t, u;
+    acb_init(t);
+    acb_init(u);
+    arb_set(acb_realref(t), nu);
+    arb_set(acb_realref(u), z);
+    acb_hypgeom_bessel_y(t, t, u, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+}
+
+void
+arb_hypgeom_bessel_jy(arb_t res1, arb_t res2, const arb_t nu, const arb_t z, slong prec)
+{
+    acb_t t, u;
+    acb_init(t);
+    acb_init(u);
+    arb_set(acb_realref(t), nu);
+    arb_set(acb_realref(u), z);
+    acb_hypgeom_bessel_jy(t, u, t, u, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res1, acb_realref(t));
+    else
+        arb_indeterminate(res1);
+    if (acb_is_finite(u) && acb_is_real(u))
+        arb_swap(res2, acb_realref(u));
+    else
+        arb_indeterminate(res2);
+    acb_clear(t);
+    acb_clear(u);
+}
+void
+arb_hypgeom_bessel_i(arb_t res, const arb_t nu, const arb_t z, slong prec)
+{
+    acb_t t, u;
+    acb_init(t);
+    acb_init(u);
+    arb_set(acb_realref(t), nu);
+    arb_set(acb_realref(u), z);
+    acb_hypgeom_bessel_i(t, t, u, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+}
+
+void
+arb_hypgeom_bessel_k(arb_t res, const arb_t nu, const arb_t z, slong prec)
+{
+    acb_t t, u;
+    acb_init(t);
+    acb_init(u);
+    arb_set(acb_realref(t), nu);
+    arb_set(acb_realref(u), z);
+    acb_hypgeom_bessel_k(t, t, u, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+}
+
diff --git a/doc/source/arb_hypgeom.rst b/doc/source/arb_hypgeom.rst
index a45849c6..8dbb2ed9 100644
--- a/doc/source/arb_hypgeom.rst
+++ b/doc/source/arb_hypgeom.rst
@@ -191,6 +191,31 @@ Exponential and trigonometric integrals
     Computes the logarithmic integral (optionally the offset version)
     of the power series *z*, truncated to length *len*.
 
+Bessel functions
+-------------------------------------------------------------------------------
+
+.. function:: void arb_hypgeom_bessel_j(arb_t res, const arb_t nu, const arb_t z, slong prec)
+
+    Computes the Bessel function of the first kind `J_{\nu}(z)`.
+
+.. function:: void arb_hypgeom_bessel_y(arb_t res, const arb_t nu, const arb_t z, slong prec)
+
+    Computes the Bessel function of the second kind `Y_{\nu}(z)`.
+
+.. function:: void arb_hypgeom_bessel_jy(arb_t res1, arb_t res2, const arb_t nu, const arb_t z, slong prec)
+
+    Sets *res1* to `J_{\nu}(z)` and *res2* to `Y_{\nu}(z)`, computed
+    simultaneously.
+
+.. function:: void arb_hypgeom_bessel_i(arb_t res, const arb_t nu, const arb_t z, slong prec)
+
+    Computes the modified Bessel function of the first kind
+    `I_{\nu}(z) = z^{\nu} (iz)^{-\nu} J_{\nu}(iz)`.
+
+.. function:: void arb_hypgeom_bessel_k(arb_t res, const arb_t nu, const arb_t z, slong prec)
+
+    Computes the modified Bessel function of the second kind `K_{\nu}(z)`.
+
 Airy functions
 -------------------------------------------------------------------------------