diff --git a/CMakeLists.txt b/CMakeLists.txt
index f178408..00b13bc 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -100,6 +100,17 @@ endif()
 include(Util)
 find_or_fetch(libneo)
 
+### fortnum: numerical core replacing the former GSL routines. libneo may pull
+### fortnum in transitively, so guard against declaring the target twice.
+include(FetchContent)
+if(NOT TARGET fortnum)
+  FetchContent_Declare(fortnum
+    GIT_REPOSITORY https://github.com/lazy-fortran/fortnum.git
+    GIT_TAG 92de6e949a772cfffc73bb5295fe5e2b056b9c18
+  )
+  FetchContent_MakeAvailable(fortnum)
+endif()
+
 ### Project source files
 # SuiteSparse
 set(SUITESPARSE_SRC_FILES
@@ -152,7 +163,7 @@ if(WITH_MFEM)
   )
 endif()
 set_source_files_properties(src/mephit_run.c ${MEPHIT_C_SRC_FILES} ${MEPHIT_CPP_SRC_FILES}
-	PROPERTIES COMPILE_FLAGS "${CMAKE_C_FLAGS} -DREAL=double -I${TRIANGLE_INCLUDE_DIR} -I${SUITESPARSE_INCLUDE_DIRS} -L${TRIANGLE_LIB_DIR}")
+	PROPERTIES COMPILE_FLAGS "${CMAKE_C_FLAGS} -DREAL=double -I${TRIANGLE_INCLUDE_DIR} -I${SUITESPARSE_INCLUDE_DIRS} -I${fortnum_SOURCE_DIR}/include -L${TRIANGLE_LIB_DIR}")
 
 
 ### Define library
@@ -176,6 +187,7 @@ set(MEPHIT_LIBS
   hdf5::hdf5_fortran
   hdf5::hdf5_hl_fortran
   GSL::gsl
+  fortnum
   BLAS::BLAS
   LAPACK::LAPACK
   ${FFTW_LIBRARIES}
diff --git a/src/mephit_pert.f90 b/src/mephit_pert.f90
index d50f372..4605bbb 100644
--- a/src/mephit_pert.f90
+++ b/src/mephit_pert.f90
@@ -1,7 +1,6 @@
 module mephit_pert
 
   use iso_fortran_env, only: dp => real64
-  use iso_c_binding, only: c_int, c_double
 
   implicit none
 
@@ -37,17 +36,6 @@ function vector_element_projection(ktri, weight, R, Z, n, f)
     end function vector_element_projection
   end interface
 
-  interface
-    function gsl_sf_bessel_icn_array(nmin, nmax, x, result_array) &
-      bind(C, name='gsl_sf_bessel_In_array')
-      import :: c_int, c_double
-      integer(c_int), intent(in), value :: nmin, nmax
-      real(c_double), intent(in), value :: x
-      real(c_double), intent(inout), dimension(nmax - nmin + 1) :: result_array
-      integer(c_int) :: gsl_sf_bessel_icn_array
-    end function gsl_sf_bessel_icn_array
-  end interface
-
   type :: L1_t
     !> Number of points on which the L1 are defined
     integer :: npoint
@@ -1427,21 +1415,21 @@ end subroutine compute_kilca_vac_coeff
   !> field
   subroutine kilca_vacuum_fourier(tor_mode, pol_mode, R_0, r, vac_coeff, &
     B_rad, B_pol, B_tor)
-    use mephit_conf, only: logger
     use mephit_util, only: imun
+    use fortnum_special, only: bessel_in_array
     integer, intent(in) :: tor_mode, pol_mode
     real(dp), intent(in) :: R_0, r
     complex(dp), intent(in) :: vac_coeff
     complex(dp), intent(out) :: B_rad, B_pol, B_tor
-    real(c_double) :: I_m(-1:1), k_z_r
-    integer(c_int) :: status
+    real(dp) :: I_n(0:abs(pol_mode) + 1), I_m(-1:1), k_z_r
 
     k_z_r = tor_mode / R_0 * r
-    status = gsl_sf_bessel_icn_array(abs(pol_mode)-1, abs(pol_mode)+1, k_z_r, I_m)
-    if (status /= 0 .and. logger%err) then
-      write (logger%msg, '("gsl_sf_bessel_In_array returned error ", i0)') status
-      call logger%write_msg
-    end if
+    ! fortnum fills I_n(0:nmax) with I_0..I_nmax; recover the order |pol_mode|-1
+    ! (which is -1 for pol_mode = 0) via the symmetry I_{-n}(x) = I_n(x).
+    call bessel_in_array(abs(pol_mode) + 1, k_z_r, I_n)
+    I_m(-1) = I_n(abs(abs(pol_mode) - 1))
+    I_m(0) = I_n(abs(pol_mode))
+    I_m(1) = I_n(abs(pol_mode) + 1)
     B_rad = 0.5d0 * (I_m(-1) + I_m(1)) * vac_coeff
     B_pol = imun * pol_mode / k_z_r * I_m(0) * vac_coeff
     B_tor = imun * I_m(0) * vac_coeff