diff --git a/sources/evolution.F90 b/sources/evolution.F90
index 5a1113b..6f02e67 100644
--- a/sources/evolution.F90
+++ b/sources/evolution.F90
@@ -1727,9 +1727,6 @@ module evolution
     integer      :: l, n, status
     real(kind=8) :: t, ds
 
-    real(kind=8), parameter :: f21 = 3.0d+00 / 4.0d+00, f22 = 1.0d+00 / 4.0d+00
-    real(kind=8), parameter :: f31 = 1.0d+00 / 3.0d+00, f32 = 2.0d+00 / 3.0d+00
-
 !-------------------------------------------------------------------------------
 !
     status = 1
@@ -1771,8 +1768,8 @@ module evolution
 
 != 2nd step: U(2) = 3/4 U(n) + 1/4 U(1) + 1/4 dt F[U(1)]
 !
-      ds = f22 * dt
-      t  = time + 0.5d+00 * dt
+      ds = dt / 4.0d+00
+      t  = time + 5.0d-01 * dt
 
 !$omp parallel do default(shared) private(pdata)
       do l = 1, n
@@ -1791,8 +1788,9 @@ module evolution
       do l = 1, n
         pdata => data_blocks(l)%ptr
 
-        pdata%uu(:,:,:,:,2) = f21 * pdata%uu(:,:,:,:,1) &
-                            + f22 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
+        pdata%uu(:,:,:,:,2) = (3.0d+00 * pdata%uu(:,:,:,:,1) &
+                                       + pdata%uu(:,:,:,:,2) &
+                                  + ds * pdata%du(:,:,:,:)) / 4.0d+00
       end do
 !$omp end parallel do
 
@@ -1802,7 +1800,7 @@ module evolution
 
 != 3rd step: U(n+1) = 1/3 U(n) + 2/3 U(2) + 2/3 dt F[U(2)]
 !
-      ds = f32 * dt
+      ds = 2.0d+00 * dt / 3.0d+00
       t  = time + dt
 
 !$omp parallel do default(shared) private(pdata)
@@ -1822,8 +1820,9 @@ module evolution
       do l = 1, n
         pdata => data_blocks(l)%ptr
 
-        pdata%uu(:,:,:,:,2) = f31 * pdata%uu(:,:,:,:,1) &
-                            + f32 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
+        pdata%uu(:,:,:,:,2) = (           pdata%uu(:,:,:,:,1) &
+                             + 2.0d+00 * (pdata%uu(:,:,:,:,2) &
+                                  + ds * pdata%du(:,:,:,:))) / 3.0d+00
 
         pdata%u => pdata%uu(:,:,:,:,1)
       end do
@@ -1845,7 +1844,7 @@ module evolution
 
       100 continue
 
-      if (status /= 0) dt = 2.5d-01 * dt
+      if (status /= 0) dt = dt / 4.0d+00
 
     end do
 
@@ -1895,9 +1894,6 @@ module evolution
     integer      :: l, n, status
     real(kind=8) :: t, ds
 
-    real(kind=8), parameter :: b1 = 2.0d+00 / 3.0d+00, b2 = 1.0d+00 / 3.0d+00, &
-                               b3 = 1.0d+00 / 6.0d+00
-
 !-------------------------------------------------------------------------------
 !
     status = 1
@@ -1988,9 +1984,9 @@ module evolution
       do l = 1, n
         pdata => data_blocks(l)%ptr
 
-        pdata%uu(:,:,:,:,2) = b1 * pdata%uu(:,:,:,:,1) &
-                            + b2 * pdata%uu(:,:,:,:,2) &
-                            + b3 * dt * pdata%du(:,:,:,:)
+        pdata%uu(:,:,:,:,2) = (2.0d+00 * pdata%uu(:,:,:,:,1) &
+                                       + pdata%uu(:,:,:,:,2) &
+                             + 5.0d-01 * dt * pdata%du(:,:,:,:)) / 3.0d+00
       end do
 !$omp end parallel do