diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,4 +1,11 @@
 
+Version 5.5
+-----
+
+* Added the 4th order Runge-Kutta 3/8-method.
+
+* Added experimental functions delayByDT and delayIByDT.
+
 Version 5.4
 -----
 
diff --git a/Simulation/Aivika/Internal/Specs.hs b/Simulation/Aivika/Internal/Specs.hs
--- a/Simulation/Aivika/Internal/Specs.hs
+++ b/Simulation/Aivika/Internal/Specs.hs
@@ -29,7 +29,8 @@
         integStopPoint,
         simulationStopPoint,
         timeGrid,
-        pointAt) where
+        pointAt,
+        delayPoint) where
 
 import Data.IORef
 
@@ -49,6 +50,7 @@
 data Method = Euler          -- ^ Euler's method
             | RungeKutta2    -- ^ the 2nd order Runge-Kutta method
             | RungeKutta4    -- ^ the 4th order Runge-Kutta method
+            | RungeKutta4b   -- ^ the 4th order Runge-Kutta 3/8-method
             deriving (Eq, Ord, Show)
 
 -- | It indentifies the simulation run.
@@ -134,6 +136,7 @@
     Euler -> [0]
     RungeKutta2 -> [0, 1]
     RungeKutta4 -> [0, 1, 2, 3]
+    RungeKutta4b -> [0, 1, 2, 3]
 
 -- | Returns the first and last integration phases.
 integPhaseBnds :: Specs -> (Int, Int)
@@ -142,6 +145,7 @@
     Euler -> (0, 0)
     RungeKutta2 -> (0, 1)
     RungeKutta4 -> (0, 3)
+    RungeKutta4b -> (0, 3)
 
 -- | Returns the first integration phase, i.e. zero.
 integPhaseLoBnd :: Specs -> Int
@@ -154,6 +158,7 @@
     Euler -> 0
     RungeKutta2 -> 1
     RungeKutta4 -> 3
+    RungeKutta4b -> 3
 
 -- | Returns a simulation time for the integration point specified by 
 -- the specs, iteration and phase.
@@ -171,6 +176,10 @@
             delta RungeKutta4 1 = spcDT sc / 2
             delta RungeKutta4 2 = spcDT sc / 2
             delta RungeKutta4 3 = spcDT sc
+            delta RungeKutta4b 0 = 0
+            delta RungeKutta4b 1 = spcDT sc / 3
+            delta RungeKutta4b 2 = 2 * spcDT sc / 3
+            delta RungeKutta4b 3 = spcDT sc
 
 -- | Return the integration time values.
 integTimes :: Specs -> [Double]
@@ -256,3 +265,23 @@
           | i == n'   = (i, t2)
           | otherwise = (i, t0 + (fromIntegral i) * dt)
   in map f [0 .. n']
+
+-- | Delay the point by the specified positive number of iterations.
+delayPoint :: Point -> Int -> Point
+delayPoint p dn
+  | dn <= 0   = error "Expected the positive number of iterations: delayPoint"
+  | otherwise =
+    let sc = pointSpecs p
+        n  = pointIteration p
+        ph = pointPhase p
+    in if ph < 0
+       then let t' = pointTime p - fromIntegral dn * spcDT sc
+                n' = fromIntegral $ floor $ (t' - spcStartTime sc) / spcDT sc
+            in p { pointTime = t',
+                   pointIteration = n',
+                   pointPhase = -1 }
+       else let n' = n - dn
+                t' = basicTime sc n' ph
+            in p { pointTime = t',
+                   pointIteration = n',
+                   pointPhase = ph }
diff --git a/Simulation/Aivika/SystemDynamics.hs b/Simulation/Aivika/SystemDynamics.hs
--- a/Simulation/Aivika/SystemDynamics.hs
+++ b/Simulation/Aivika/SystemDynamics.hs
@@ -49,6 +49,8 @@
         -- * Discrete Functions
         delay,
         delayI,
+        delayByDT,
+        delayIByDT,
         step,
         pulse,
         pulseP,
@@ -237,6 +239,79 @@
     _ -> 
       error "Incorrect phase: integRK4"
 
+integRK4b :: Dynamics Double
+             -> Dynamics Double
+             -> Dynamics Double
+             -> Point -> IO Double
+integRK4b (Dynamics f) (Dynamics i) (Dynamics y) p =
+  case pointPhase p of
+    0 -> case pointIteration p of
+      0 -> 
+        i p
+      n -> do
+        let sc = pointSpecs p
+            ty = basicTime sc (n - 1) 0
+            t1 = ty
+            t2 = basicTime sc (n - 1) 1
+            t3 = basicTime sc (n - 1) 2
+            t4 = basicTime sc (n - 1) 3
+            py = p { pointTime = ty, pointIteration = n - 1, pointPhase = 0 }
+            p1 = py
+            p2 = p { pointTime = t2, pointIteration = n - 1, pointPhase = 1 }
+            p3 = p { pointTime = t3, pointIteration = n - 1, pointPhase = 2 }
+            p4 = p { pointTime = t4, pointIteration = n - 1, pointPhase = 3 }
+        vy <- y py
+        k1 <- f p1
+        k2 <- f p2
+        k3 <- f p3
+        k4 <- f p4
+        let !v = vy + spcDT sc / 8.0 * (k1 + 3.0 * (k2 + k3) + k4)
+        return v
+    1 -> do
+      let sc = pointSpecs p
+          n  = pointIteration p
+          ty = basicTime sc n 0
+          t1 = ty
+          py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }
+          p1 = py
+      vy <- y py
+      k1 <- f p1
+      let !v = vy + spcDT sc / 3.0 * k1
+      return v
+    2 -> do
+      let sc = pointSpecs p
+          n  = pointIteration p
+          ty = basicTime sc n 0
+          t1 = ty
+          t2 = basicTime sc n 1
+          py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }
+          p1 = py
+          p2 = p { pointTime = t2, pointIteration = n, pointPhase = 1 }
+      vy <- y py
+      k1 <- f p1
+      k2 <- f p2
+      let !v = vy + spcDT sc * (- k1 / 3.0 + k2)
+      return v
+    3 -> do
+      let sc = pointSpecs p
+          n  = pointIteration p
+          ty = basicTime sc n 0
+          t1 = ty
+          t2 = basicTime sc n 1
+          t3 = basicTime sc n 2
+          py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }
+          p1 = py
+          p2 = p { pointTime = t2, pointIteration = n, pointPhase = 1 }
+          p3 = p { pointTime = t3, pointIteration = n, pointPhase = 2 }
+      vy <- y py
+      k1 <- f p1
+      k2 <- f p2
+      k3 <- f p3
+      let !v = vy + spcDT sc * (k1 - k2 + k3)
+      return v
+    _ -> 
+      error "Incorrect phase: integRK4b"
+
 -- | Return an integral with the specified derivative and initial value.
 --
 -- To create a loopback, you should use the recursive do-notation.
@@ -263,6 +338,7 @@
           Euler -> return $ Dynamics $ integEuler diff i y
           RungeKutta2 -> return $ Dynamics $ integRK2 diff i y
           RungeKutta4 -> return $ Dynamics $ integRK4 diff i y
+          RungeKutta4b -> return $ Dynamics $ integRK4b diff i y
       return y
 
 integEulerEither :: Dynamics (Either Double Double)
@@ -632,10 +708,71 @@
                                   pointIteration = n',
                                   pointPhase = -1 }
             | n' > n    = error $
-                          "Cannot return the future data: delay. " ++
+                          "Cannot return the future data: delayI. " ++
                           "The lag time cannot be negative."
             | otherwise = error $
-                          "Cannot return the current data: delay. " ++
+                          "Cannot return the current data: delayI. " ++
+                          "The lag time is too small."
+      y
+
+-- | Return the delayed value by the specified positive number of
+-- integration time steps used for calculating the lag time.
+delayByDT :: Dynamics a
+             -- ^ the value to delay
+             -> Dynamics Int
+             -- ^ the delay as a multiplication of the corresponding number
+             -- and the integration time step
+             -> Dynamics a
+             -- ^ the delayed value
+delayByDT (Dynamics x) (Dynamics d) = discreteDynamics $ Dynamics r 
+  where
+    r p = do 
+      let sc = pointSpecs p
+          n  = pointIteration p
+      a <- d p
+      let p' = delayPoint p a
+          n' = pointIteration p'
+          y | n' < 0    = x $ p { pointTime = spcStartTime sc,
+                                  pointIteration = 0, 
+                                  pointPhase = 0 }
+            | n' < n    = x p'
+            | n' > n    = error $
+                          "Cannot return the future data: delayByDT. " ++
+                          "The lag time cannot be negative."
+            | otherwise = error $
+                          "Cannot return the current data: delayByDT. " ++
+                          "The lag time is too small."
+      y
+      
+-- | Return the delayed value by the specified initial value and
+-- a positive number of integration time steps used for calculating
+-- the lag time. It allows creating a loop back.
+delayIByDT :: Dynamics a
+              -- ^ the value to delay
+              -> Dynamics Int
+              -- ^ the delay as a multiplication of the corresponding number
+              -- and the integration time step
+              -> Dynamics a
+              -- ^ the initial value
+              -> Simulation (Dynamics a)
+              -- ^ the delayed value
+delayIByDT (Dynamics x) (Dynamics d) (Dynamics i) = M.memoDynamics $ Dynamics r 
+  where
+    r p = do 
+      let sc = pointSpecs p
+          n  = pointIteration p
+      a <- d p
+      let p' = delayPoint p a
+          n' = pointIteration p'
+          y | n' < 0    = i $ p { pointTime = spcStartTime sc,
+                                  pointIteration = 0, 
+                                  pointPhase = 0 }
+            | n' < n    = x p'
+            | n' > n    = error $
+                          "Cannot return the future data: delayIByDT. " ++
+                          "The lag time cannot be negative."
+            | otherwise = error $
+                          "Cannot return the current data: delayIByDT. " ++
                           "The lag time is too small."
       y
 
diff --git a/aivika.cabal b/aivika.cabal
--- a/aivika.cabal
+++ b/aivika.cabal
@@ -1,5 +1,5 @@
 name:            aivika
-version:         5.4
+version:         5.5
 synopsis:        A multi-method simulation library
 description:
     Aivika is a discrete event simulation (DES) framework with support of activity-oriented,