diff --git a/sources/evolution.F90 b/sources/evolution.F90
index 0a55fc2..8e417ab 100644
--- a/sources/evolution.F90
+++ b/sources/evolution.F90
@@ -83,6 +83,9 @@ module evolution
 !
   real(kind=8)      , save :: atol        = 1.0d-04
   real(kind=8)      , save :: rtol        = 1.0d-04
+  real(kind=8)      , save :: fac         = 9.0d-01
+  real(kind=8)      , save :: facmin      = 1.0d-01
+  real(kind=8)      , save :: facmax      = 5.0d+00
   integer           , save :: mrej        = 5
   integer           , save :: niterations = 0
   integer           , save :: nrejections = 0
@@ -190,11 +193,14 @@ module evolution
     call get_parameter("glm_alpha"   , glm_alpha  )
 
 ! get the absolute and relative tolerances, the maximum number of passes for
-! the adaptive step
+! the adaptive step, and limiting factors
 !
     call get_parameter("absolute_tolerance", atol)
     call get_parameter("relative_tolerance", rtol)
     call get_parameter("maximum_rejections", mrej)
+    call get_parameter("factor"            , fac )
+    call get_parameter("minimum_factor"    , facmin)
+    call get_parameter("maximum_factor"    , facmax)
 
 ! select the integration method, check the correctness of the integration
 ! parameters and adjust the CFL coefficient if necessary
@@ -426,6 +432,9 @@ module evolution
       call print_parameter(verbose, "absolute tolerance", atol)
       call print_parameter(verbose, "relative tolerance", rtol)
       call print_parameter(verbose, "maximum rejections", mrej)
+      call print_parameter(verbose, "factor"            , fac)
+      call print_parameter(verbose, "minimum factor"    , facmin)
+      call print_parameter(verbose, "maximum factor"    , facmax)
 
     end if
 
@@ -954,11 +963,14 @@ module evolution
     logical      :: test
     integer      :: n, l, nrej
     real(kind=8) :: tm, dtm, ds, umax, utol, emax
+    real(kind=8) :: fc, fcmn, fcmx
 
     logical     , save :: first = .true.
     real(kind=8), save :: ft, fl, fr, gt, gl
 
     real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr
+
+    real(kind=8), parameter :: k1 = -2.9d-01, k2 = 1.05d-01, k3 = -5.0d-02
 !
 !-------------------------------------------------------------------------------
 !
@@ -979,6 +991,10 @@ module evolution
 
     end if
 
+    fc   = fac
+    fcmn = facmin
+    fcmx = facmax
+
     l = get_dblocks()
 #if NDIMS == 3
     allocate(lerr(l,nf,nc,nc,nc))
@@ -1105,15 +1121,19 @@ module evolution
         errs(2) = errs(1)
         errs(1) = emax
 
-        dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00)                 &
-                                       * errs(3)**(-0.025d+00)
+        dte = dt * min(fcmx, max(fcmn,                                         &
+                                 fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))
+
+        fcmx = facmax
       else
         errs(1) = emax
 
-        dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00)                 &
-                                       * errs(3)**(-0.025d+00)
+        dte = dt * min(fcmx, max(fcmn,                                         &
+                                 fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))
         dt  = dte
 
+        fcmx = fac
+
         nrej        = nrej        + 1 ! rejection count in the current step
         nrejections = nrejections + 1
       end if
@@ -1944,8 +1964,12 @@ module evolution
     logical      :: test
     integer      :: i, l, nrej
     real(kind=8) :: tm, dtm, ds, umax, utol, emax
+    real(kind=8) :: fc, fcmn, fcmx
 
     real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr
+
+    real(kind=8), parameter :: k1 = -5.8d-01 / 3.0d+00, k2 = 7.0d-02,          &
+                               k3 = -1.0d-01 / 3.0d+00
 !
 !-------------------------------------------------------------------------------
 !
@@ -1973,6 +1997,10 @@ module evolution
 
     end if
 
+    fc   = fac
+    fcmn = facmin
+    fcmx = facmax
+
     l = get_dblocks()
 #if NDIMS == 3
     allocate(lerr(l,nf,nc,nc,nc))
@@ -2179,15 +2207,19 @@ module evolution
         errs(2) = errs(1)
         errs(1) = emax
 
-        dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00)                 &
-                                       * errs(3)**(-0.025d+00)
+        dte = dt * min(fcmx, max(fcmn,                                         &
+                                 fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))
+
+        fcmx = facmax
       else
         errs(1) = emax
 
-        dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00)                 &
-                                       * errs(3)**(-0.025d+00)
+        dte = dt * min(fcmx, max(fcmn,                                         &
+                                 fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))
         dt  = dte
 
+        fcmx = fac
+
         nrej        = nrej        + 1 ! rejection count in the current step
         nrejections = nrejections + 1
       end if