diff --git a/input/gray.ini b/input/gray.ini
index 60a0aee..20929c4 100644
--- a/input/gray.ini
+++ b/input/gray.ini
@@ -26,6 +26,14 @@ nstep = 12000
 ;  2: real part of the eikonal (S_r)
 idst = 0
 
+; Choice of the integration method
+;   0: Explicit Euler (1⁰ order)
+;   1: Semi-implicit Euler (1⁰ order, symplectic)
+;   2: Velocity Verlet (2⁰ order, symplectic)
+;   3: 2-stage Runge-Kutta (2⁰ order)
+;   4: 4-stage Runge-Kutta (4⁰ order)
+integrator = 4
+
 
 [ecrh_cd]
 ; Choice of the power absorption model
diff --git a/src/beamdata.f90 b/src/beamdata.f90
index ee8c768..a36c84a 100644
--- a/src/beamdata.f90
+++ b/src/beamdata.f90
@@ -3,14 +3,14 @@ module beamdata
   implicit none
 
   integer, save :: nray,nrayr,nrayth,nstep,jkray1
-  real(wp_), save :: dst,h,hh,h6,rwmax,twodr2
+  real(wp_), save :: dst,rwmax,twodr2
 
 contains
 
   subroutine init_btr(rtrparam,ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
                       tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
     use gray_params, only : raytracing_parameters
-    use const_and_precisions, only : half,two
+    use const_and_precisions, only : two
     implicit none
     type(raytracing_parameters), intent(inout) :: rtrparam
     real(wp_),  dimension(:,:), intent(out), pointer :: ywork,ypwork, &
@@ -24,10 +24,6 @@ contains
     integer :: jray1
 
     dst = rtrparam%dst
-    h   = dst
-    hh  = h*half
-    h6  = h/6.0_wp_
-    
     nrayr  = rtrparam%nrayr
     nrayth = rtrparam%nrayth
     rwmax  = rtrparam%rwmax
diff --git a/src/gray_core.f90 b/src/gray_core.f90
index 5f3468b..d8a1040 100644
--- a/src/gray_core.f90
+++ b/src/gray_core.f90
@@ -5,6 +5,15 @@ module gray_core
 
   implicit none
 
+  abstract interface
+    function rhs_function(y) result(f)
+      ! Function passed to the integrator subroutine
+      use const_and_precisions, only : wp_
+      real(wp_), intent(in) :: y(6)
+      real(wp_)             :: f(6)
+    end function
+  end interface
+
 contains
 
   subroutine gray_main(params, data, results, error, rhout)
@@ -265,7 +274,9 @@ contains
           do jk=1,params%raytracing%nray
             if(iwait(jk)) cycle                                                            ! jk ray is waiting for next pass
             stv(jk) = stv(jk) + params%raytracing%dst                                      ! current ray step
-            call rkstep(sox,bres,xgcn,yw(:,jk),ypw(:,jk),gri(:,jk),ggri(:,:,jk),igrad_b)
+            call integrator(yw(:, jk), ypw(:, jk), rhs, &
+                            params%raytracing%dst,      &
+                            params%raytracing%integrator)
           end do
           ! update position and grad
           if(igrad_b == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri)
@@ -613,6 +624,39 @@ contains
        stnext,stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam,     &
        pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv)
     ! =========== free memory END ===========
+
+    contains
+
+      ! Functions that needs the scope of gray_main
+
+      function rhs(y) result(dery)
+        ! Computes the right-hand-side terms of the ray equations
+        ! To be passed to the integrator subroutine
+        implicit none
+
+        ! function arguments
+        real(wp_), dimension(6), intent(in) :: y ! (x̅, N̅)
+        ! result
+        real(wp_), dimension(6) :: dery
+
+        ! local variables
+        real(wp_) :: xg, yg
+        real(wp_), dimension(3) :: xv, anv, bv, derxg, deryg
+        real(wp_), dimension(3, 3) :: derbv
+
+        xv  = y(1:3)  ! x̅ position
+        anv = y(4:6)  ! N̅ refractive index
+
+        ! computes derivatives of plasma quantities: B̅, ∇X, ∇Y, ∇B̅
+        ! (these are needed for the next part)
+        call plas_deriv(xv, bres, xgcn, dens, btot, bv, derbv, &
+                        xg, yg, derxg, deryg)
+
+        ! computes derivatives of dispersion relation: ∂Λ/∂x̅, ∂Λ/∂N̅
+        call disp_deriv(anv, sox, xg, yg, derxg, deryg, bv, derbv, &
+                        gri(:, jk), ggri(:, :, jk), igrad_b, dery)
+      end function rhs
+
   end subroutine gray_main
 
 
@@ -966,55 +1010,65 @@ contains
 
 
 
-  subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr,igrad)
-!     Runge-Kutta integrator
-!    use gray_params, only : igrad
-    use beamdata, only : h,hh,h6
+  subroutine integrator(y, yp, f, h, method)
+    ! Integrator of the raytracing equations
     implicit none
-    real(wp_), intent(in) :: bres,xgcn
-    real(wp_), dimension(6), intent(inout) :: y
-    real(wp_), dimension(6), intent(in) :: yp
-    real(wp_), dimension(3), intent(in) :: dgr
-    real(wp_), dimension(3,3), intent(in) :: ddgr
-    integer, intent(in) :: igrad,sox
 
-    real(wp_), dimension(6) :: yy,fk1,fk2,fk3,fk4
+    real(wp_), dimension(6), intent(inout) :: y  ! y̅  = (x̅, N̅)
+    real(wp_), dimension(6), intent(in)    :: yp ! y̅˙ = (dx̅/dσ, dN̅/dσ)
+    procedure(rhs_function)                :: f  ! dy̅/dσ = f̅(y̅)
+    real(wp_),               intent(in)    :: h  ! step size
+    integer,                 intent(in)    :: method
 
-    fk1 = yp
-    yy = y + fk1*hh
-    call rhs(sox,bres,xgcn,yy,dgr,ddgr,fk2,igrad)
-    yy = y + fk2*hh
-    call rhs(sox,bres,xgcn,yy,dgr,ddgr,fk3,igrad)
-    yy = y + fk3*h
-    call rhs(sox,bres,xgcn,yy,dgr,ddgr,fk4,igrad)
+    ! local variables
+    real(wp_), dimension(6) :: yy, k1, k2, k3, k4
 
-    y = y + h6*(fk1 + 2*fk2 + 2*fk3 + fk4)
-  end subroutine rkstep
+    select case (method)
+      case (0)
+        ! Explicit Euler (1⁰ order)
+        y = y + yp*h
 
+      case (1)
+        ! Semi-implicit Euler (1⁰ order, symplectic)
+        ! P = p - ∂H/∂q(q, p)
+        ! Q = q + ∂H/∂p(q, P)
+        k1 = yp
+        y(4:6) = y(4:6) + k1(4:6)*h
+        k2 = f(y)
+        y(1:3) = y(1:3) + k2(1:3)*h
 
+      case (2)
+        ! Velocity Verlet (2⁰ order, symplectic)
+        ! p(n+½) = p(n) - ∂H/∂q(q(n), p(n))⋅h/2
+        ! q(n+1) = q(n) + ∂H/∂p(q(n), p(n+½))⋅h
+        ! p(n+1) = p(n+½) - ∂H/∂q(q(n+1), p(n+½))⋅h/2
+        k1 = yp
+        y(4:6) = y(4:6) + k1(4:6)*h/2
+        k2 = f(y)
+        y(1:3) = y(1:3) + k2(1:3)*h
+        k3 = f(y)
+        y(4:6) = y(4:6) + k3(4:6)*h/2
 
-  subroutine rhs(sox,bres,xgcn,y,dgr,ddgr,dery,igrad)
-! Compute right-hand side terms of the ray equations (dery)
-! used in R-K integrator
-    implicit none
-! arguments
-    real(wp_), dimension(6), intent(in) :: y
-    real(wp_), intent(in) :: bres,xgcn
-    real(wp_), dimension(3), intent(in) :: dgr
-    real(wp_), dimension(3,3), intent(in) :: ddgr
-    real(wp_), dimension(6), intent(out) :: dery
-    integer, intent(in) :: igrad,sox
-! local variables
-    real(wp_) :: dens,btot,xg,yg
-    real(wp_), dimension(3) :: xv,anv,bv,derxg,deryg
-    real(wp_), dimension(3,3) :: derbv
-
-    xv = y(1:3)
-    call plas_deriv(xv,bres,xgcn,dens,btot,bv,derbv,xg,yg,derxg,deryg)
-    anv = y(4:6)
-    call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,dgr,ddgr,igrad,dery)
-  end subroutine rhs
+      case (3)
+        ! 2-stage Runge-Kutta (2⁰ order)
+        k1 = yp
+        yy = y + k1*h
+        k2 = f(yy)
+        yy = y + k2*h
+        y  = y + (k1 + k2)*h/2
 
+      case default
+        ! 4-stage Runge-Kutta (4⁰ order)
+        k1 = yp
+        yy = y + k1*h/2
+        k2 = f(yy)
+        yy = y + k2*h/2
+        k3 = f(yy)
+        yy = y + k3*h
+        k4 = f(yy)
+        y  = y + (k1 + 2*k2 + 2*k3 + k4)*h/6
+    end select
+  end subroutine integrator
 
 
   subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, &
diff --git a/src/gray_params.f90 b/src/gray_params.f90
index 416ca85..790d483 100644
--- a/src/gray_params.f90
+++ b/src/gray_params.f90
@@ -5,6 +5,9 @@ module gray_params
   integer, parameter :: lenfnm = 256
   integer, parameter :: headw = 132, headl = 21
 
+  ! Note: when adding/removing a parameter remember to keep
+  ! the gray_params.sh script in sync with this file
+
   ! Antenna/wave launcher parameters
   type antenna_parameters
     ! From gray_params.data:
@@ -55,6 +58,7 @@ module gray_params
     integer :: igrad                ! Complex eikonal switch
     integer :: nstep                ! Max number of integration steps
     integer :: idst                 ! Choice of the integration variable
+    integer :: integrator           ! Choice of the integration method
     integer :: ipass                ! Number of plasma passes
     integer :: ipol                 ! Whether to compute wave polarisation
   end type
@@ -365,6 +369,11 @@ contains
     read(u, *) params%output%istpr, params%output%istpl
 
     close(u)
+
+    ! Default values of parameters introduced after
+    ! gray_params.data has been deprecated
+    params%raytracing%integrator = 4
+
   end subroutine read_gray_params
 
 
diff --git a/src/gray_params.sh b/src/gray_params.sh
index ae21d72..3351119 100644
--- a/src/gray_params.sh
+++ b/src/gray_params.sh
@@ -11,7 +11,7 @@ sets='antenna equilibrium profiles raytracing ecrh_cd output misc'
 antenna='alpha beta power psi chi iox ibeam filenm fghz pos w ri phi'
 equilibrium='ssplps ssplf factb sgnb sgni ixp iequil icocos ipsinorm idesc ifreefmt filenm'
 profiles='psnbnd sspld factne factte iscal irho iprof filenm'
-raytracing='rwmax dst nrayr nrayth nstep igrad idst ipass ipol'
+raytracing='rwmax dst nrayr nrayth nstep igrad idst ipass ipol integrator'
 ecrh_cd='iwarm ilarm imx ieccd'
 output='ipec nrho istpr istpl'
 misc='rwall'