diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright (c) 2009, 2010, 2011 David Sorokin <david.sorokin@gmail.com>
+Copyright (c) 2009, 2010, 2011, 2012 David Sorokin <david.sorokin@gmail.com>
 
 All rights reserved.
 
diff --git a/Simulation/Aivika/Dynamics.hs b/Simulation/Aivika/Dynamics.hs
--- a/Simulation/Aivika/Dynamics.hs
+++ b/Simulation/Aivika/Dynamics.hs
@@ -9,27 +9,12 @@
 --
 -- The module defines the 'Dynamics' monad representing an abstract dynamic 
 -- process, i.e. a time varying polymorphic function. 
--- 
--- This is a key point of the Aivika simulation library. With help of this monad 
--- we can simulate the system of ordinary differential equations (ODEs) of 
--- System Dynamics, define the tasks of Discrete Event Simulation (DES) supporting 
--- different paradigms. Also we can use the Agent-based Modeling. Thus, 
--- we can create hybrid simulation models.
 --
 module Simulation.Aivika.Dynamics 
        (Dynamics,
-        Specs(..),
-        Method(..),
-        runDynamics1,
-        runDynamics1_,
-        runDynamics,
-        runDynamics_,
-        runDynamicsIO,
-        runDynamicsSeries1,
-        runDynamicsSeries1_,
-        runDynamicsSeries,
-        runDynamicsSeries_,
-        printDynamics1,
-        printDynamics) where
+        DynamicsLift(..),
+        runDynamicsInStart,
+        runDynamicsInFinal,
+        runDynamics) where
 
 import Simulation.Aivika.Dynamics.Internal.Dynamics
diff --git a/Simulation/Aivika/Dynamics/Agent.hs b/Simulation/Aivika/Dynamics/Agent.hs
--- a/Simulation/Aivika/Dynamics/Agent.hs
+++ b/Simulation/Aivika/Dynamics/Agent.hs
@@ -30,6 +30,7 @@
 import Data.IORef
 import Control.Monad
 
+import Simulation.Aivika.Dynamics.Internal.Simulation
 import Simulation.Aivika.Dynamics.Internal.Dynamics
 import Simulation.Aivika.Dynamics.EventQueue
 
@@ -116,11 +117,13 @@
 addTimeout :: AgentState -> Double -> Dynamics () -> Dynamics ()
 addTimeout st dt (Dynamics action) =
   Dynamics $ \p ->
-  do v <- readIORef (stateVersionRef st)
+  do let q = agentQueue (stateAgent st)
+         Dynamics m0 = queueRun q
+     m0 p    -- ensure that the agent state is actual
+     v <- readIORef (stateVersionRef st)
      let m1 = Dynamics $ \p ->
            do v' <- readIORef (stateVersionRef st)
               when (v == v') $ action p
-         q = agentQueue (stateAgent st)
          Dynamics m2 = enqueue q (pointTime p + dt) m1
      m2 p
 
@@ -130,11 +133,13 @@
 addTimer :: AgentState -> Dynamics Double -> Dynamics () -> Dynamics ()
 addTimer st (Dynamics dt) (Dynamics action) =
   Dynamics $ \p ->
-  do v <- readIORef (stateVersionRef st)
+  do let q = agentQueue (stateAgent st)
+         Dynamics m0 = queueRun q
+     m0 p    -- ensure that the agent state is actual
+     v <- readIORef (stateVersionRef st)
      let m1 = Dynamics $ \p ->
            do v' <- readIORef (stateVersionRef st)
               when (v == v') $ do { m2 p; action p }
-         q = agentQueue (stateAgent st)
          Dynamics m2 = 
            Dynamics $ \p ->
            do dt' <- dt p
@@ -143,9 +148,9 @@
      m2 p
 
 -- | Create a new state.
-newState :: Agent -> Dynamics AgentState
+newState :: Agent -> Simulation AgentState
 newState agent =
-  Dynamics $ \p ->
+  Simulation $ \r ->
   do aref <- newIORef $ return ()
      dref <- newIORef $ return ()
      vref <- newIORef 0
@@ -156,9 +161,9 @@
                          stateVersionRef = vref }
 
 -- | Create a child state.
-newSubstate :: AgentState -> Dynamics AgentState
+newSubstate :: AgentState -> Simulation AgentState
 newSubstate parent =
-  Dynamics $ \p ->
+  Simulation $ \r ->
   do let agent = stateAgent parent 
      aref <- newIORef $ return ()
      dref <- newIORef $ return ()
@@ -170,9 +175,9 @@
                          stateVersionRef = vref }
 
 -- | Create an agent bound with the specified event queue.
-newAgent :: EventQueue -> Dynamics Agent
+newAgent :: EventQueue -> Simulation Agent
 newAgent queue =
-  Dynamics $ \p ->
+  Simulation $ \r ->
   do modeRef    <- newIORef CreationMode
      stateRef   <- newIORef Nothing
      return Agent { agentQueue = queue,
@@ -247,14 +252,14 @@
          "the state activation: initState."
 
 -- | Set the activation computation for the specified state.
-stateActivation :: AgentState -> Dynamics () -> Dynamics ()
+stateActivation :: AgentState -> Dynamics () -> Simulation ()
 stateActivation st action =
-  Dynamics $ \p ->
+  Simulation $ \r ->
   writeIORef (stateActivateRef st) action
   
 -- | Set the deactivation computation for the specified state.
-stateDeactivation :: AgentState -> Dynamics () -> Dynamics ()
+stateDeactivation :: AgentState -> Dynamics () -> Simulation ()
 stateDeactivation st action =
-  Dynamics $ \p ->
+  Simulation $ \r ->
   writeIORef (stateDeactivateRef st) action
   
diff --git a/Simulation/Aivika/Dynamics/Base.hs b/Simulation/Aivika/Dynamics/Base.hs
--- a/Simulation/Aivika/Dynamics/Base.hs
+++ b/Simulation/Aivika/Dynamics/Base.hs
@@ -12,12 +12,12 @@
 
 module Simulation.Aivika.Dynamics.Base
        (-- * Time Parameters
-        starttime,
+        starttime, 
         stoptime,
         dt,
         time,
         -- * Interpolation and Initial Value
-        initD,
+        initDynamics,
         discrete,
         interpolate,
         -- * Memoization
@@ -26,12 +26,12 @@
         memo0,
         umemo0,
         -- * Iterating
-        iterateD,
+        iterateDynamics,
         -- * Fold
-        foldD1,
-        foldD,
+        foldDynamics1,
+        foldDynamics,
         -- * Norming
-        divideD) where
+        divideDynamics) where
 
 import Simulation.Aivika.Dynamics.Internal.Dynamics
 import Simulation.Aivika.Dynamics.Internal.Time
diff --git a/Simulation/Aivika/Dynamics/Cont.hs b/Simulation/Aivika/Dynamics/Cont.hs
--- a/Simulation/Aivika/Dynamics/Cont.hs
+++ b/Simulation/Aivika/Dynamics/Cont.hs
@@ -7,10 +7,8 @@
 -- Stability  : experimental
 -- Tested with: GHC 7.0.3
 --
--- The 'Cont' monad looks somewhere like the standard ContT monad transformer 
--- parameterized by the 'Dynamics' monad, although this analogy is not strong. 
--- The main idea is to represent the continuation as a dynamic process varying 
--- in time.
+-- The 'Cont' monad is a variation of the standard Cont monad, where
+-- the result of applying the continuation is a dynamic process.
 --
 module Simulation.Aivika.Dynamics.Cont
        (Cont,
diff --git a/Simulation/Aivika/Dynamics/EventQueue.hs b/Simulation/Aivika/Dynamics/EventQueue.hs
--- a/Simulation/Aivika/Dynamics/EventQueue.hs
+++ b/Simulation/Aivika/Dynamics/EventQueue.hs
@@ -22,21 +22,22 @@
 import Data.IORef
 import Control.Monad
 
+import Simulation.Aivika.Dynamics.Internal.Simulation
 import Simulation.Aivika.Dynamics.Internal.Dynamics
 import qualified Simulation.Aivika.PriorityQueue as PQ
 
 -- | The 'EventQueue' type represents the event queue.
 data EventQueue = EventQueue { 
-  queuePQ   :: PQ.PriorityQueue (Dynamics (() -> IO ())),
+  queuePQ   :: PQ.PriorityQueue (() -> Dynamics ()),
   queueRun  :: Dynamics (),   -- ^ Run the event queue processing its events
   queueBusy :: IORef Bool,
   queueTime :: IORef Double }
 
 -- | Create a new event queue.
-newQueue :: Dynamics EventQueue
+newQueue :: Simulation EventQueue
 newQueue = 
-  Dynamics $ \p ->
-  do let sc = pointSpecs p
+  Simulation $ \r ->
+  do let sc = runSpecs r
      f <- newIORef False
      t <- newIORef $ spcStartTime sc
      pq <- PQ.newQueue
@@ -47,14 +48,13 @@
      return q
              
 -- | Enqueue the event which must be actuated at the specified time.
-enqueueCont :: EventQueue -> Double -> Dynamics (() -> IO ()) -> Dynamics ()
+enqueueCont :: EventQueue -> Double -> (() -> Dynamics ()) -> Dynamics ()
 enqueueCont q t c = Dynamics r where
   r p = let pq = queuePQ q in PQ.enqueue pq t c
     
 -- | Enqueue the event which must be actuated at the specified time.
 enqueue :: EventQueue -> Double -> Dynamics () -> Dynamics ()
-enqueue q t (Dynamics m) = enqueueCont q t (Dynamics c) where
-  c p = let f () = m p in return f
+enqueue q t m = enqueueCont q t (const m) 
     
 -- | Run the event queue processing its events.
 runQueue :: EventQueue -> Dynamics ()
@@ -70,7 +70,7 @@
     do let pq = queuePQ q
        f <- PQ.queueNull pq
        unless f $
-         do (t2, Dynamics c2) <- PQ.queueFront pq
+         do (t2, c2) <- PQ.queueFront pq
             let t = queueTime q
             t' <- readIORef t
             when (t2 < t') $ 
@@ -82,8 +82,8 @@
                      t0  = spcStartTime sc
                      dt  = spcDT sc
                      n2  = fromInteger $ toInteger $ floor ((t2 - t0) / dt)
-                 k <- c2 $ p { pointTime = t2,
-                               pointIteration = n2,
-                               pointPhase = -1 }
-                 k ()    -- raise the event
+                     Dynamics k = c2 ()
+                 k $ p { pointTime = t2,
+                         pointIteration = n2,
+                         pointPhase = -1 }
                  call q p
diff --git a/Simulation/Aivika/Dynamics/Internal/Cont.hs b/Simulation/Aivika/Dynamics/Internal/Cont.hs
--- a/Simulation/Aivika/Dynamics/Internal/Cont.hs
+++ b/Simulation/Aivika/Dynamics/Internal/Cont.hs
@@ -7,10 +7,8 @@
 -- Stability  : experimental
 -- Tested with: GHC 7.0.3
 --
--- The 'Cont' monad looks somewhere like the standard ContT monad transformer 
--- parameterized by the 'Dynamics' monad, although this analogy is not strong. 
--- The main idea is to represent the continuation as a dynamic process varying 
--- in time.
+-- The 'Cont' monad is a variation of the standard Cont monad, where
+-- the result of applying the continuation is a dynamic process.
 --
 module Simulation.Aivika.Dynamics.Internal.Cont
        (Cont(..),
@@ -19,20 +17,23 @@
 import Control.Monad
 import Control.Monad.Trans
 
+import Simulation.Aivika.Dynamics.Internal.Simulation
 import Simulation.Aivika.Dynamics.Internal.Dynamics
-import Simulation.Aivika.Dynamics.Lift
 
 -- | The 'Cont' type is similar to the standard Cont monad but only
--- the continuation is represented as a dynamic process varying in time.
-newtype Cont a = Cont (Dynamics (a -> IO ()) -> Dynamics ())
+-- the continuation uses a dynamic process as a result.
+newtype Cont a = Cont ((a -> Dynamics ()) -> Dynamics ())
 
 instance Monad Cont where
   return  = returnC
   m >>= k = bindC m k
 
-instance Lift Cont where
-  liftD = liftC
+instance SimulationLift Cont where
+  liftSimulation = liftSC
 
+instance DynamicsLift Cont where
+  liftDynamics = liftDC
+
 instance Functor Cont where
   fmap = liftM
 
@@ -41,44 +42,44 @@
 
 returnC :: a -> Cont a
 {-# INLINE returnC #-}
-returnC a = 
-  Cont $ \(Dynamics c) -> 
-  Dynamics $ \p -> 
-  do cont' <- c p
-     cont' a
+returnC a = Cont $ \c -> c a
                           
 bindC :: Cont a -> (a -> Cont b) -> Cont b
 {-# INLINE bindC #-}
-bindC (Cont m) k =
-  Cont $ \c ->
-  m $ Dynamics $ \p -> 
-  let cont' a = let (Cont m') = k a
-                    (Dynamics u) = m' c
-                in u p
-  in return cont'
+bindC (Cont m) k = Cont $ \c -> m (\a -> let Cont m' = k a in m' c)
 
 -- | Run the 'Cont' computation.
-runCont :: Cont a -> IO (a -> IO ()) -> Dynamics ()
+runCont :: Cont a -> (a -> Dynamics ()) -> Dynamics ()
 {-# INLINE runCont #-}
-runCont (Cont m) f = m $ Dynamics $ const f
+runCont (Cont m) = m
 
+-- | Lift the 'Simulation' computation.
+liftSC :: Simulation a -> Cont a
+{-# INLINE liftSC #-}
+liftSC (Simulation m) = 
+  Cont $ \c ->
+  Dynamics $ \p ->
+  do a <- m $ pointRun p
+     let Dynamics m' = c a
+     m' p
+     
 -- | Lift the 'Dynamics' computation.
-liftC :: Dynamics a -> Cont a
-{-# INLINE liftC #-}
-liftC (Dynamics m) =
-  Cont $ \(Dynamics c) ->
+liftDC :: Dynamics a -> Cont a
+{-# INLINE liftDC #-}
+liftDC (Dynamics m) =
+  Cont $ \c ->
   Dynamics $ \p ->
-  do cont' <- c p
-     a <- m p
-     cont' a
+  do a <- m p
+     let Dynamics m' = c a
+     m' p
      
 -- | Lift the IO computation.
 liftIOC :: IO a -> Cont a
 {-# INLINE liftIOC #-}
 liftIOC m =
-  Cont $ \(Dynamics c) ->
+  Cont $ \c ->
   Dynamics $ \p ->
-  do cont' <- c p
-     a <- m
-     cont' a
+  do a <- m
+     let Dynamics m' = c a
+     m' p
   
diff --git a/Simulation/Aivika/Dynamics/Internal/Dynamics.hs b/Simulation/Aivika/Dynamics/Internal/Dynamics.hs
--- a/Simulation/Aivika/Dynamics/Internal/Dynamics.hs
+++ b/Simulation/Aivika/Dynamics/Internal/Dynamics.hs
@@ -9,31 +9,15 @@
 --
 -- The module defines the 'Dynamics' monad representing an abstract dynamic 
 -- process, i.e. a time varying polymorphic function. 
--- 
--- This is a key point of the Aivika simulation library. With help of this monad 
--- we can simulate the system of ordinary differential equations (ODEs) of 
--- System Dynamics, define the tasks of Discrete Event Simulation (DES) supporting 
--- different paradigms. Also we can use the Agent-based Modeling. Thus, 
--- we can create hybrid simulation models.
 --
 module Simulation.Aivika.Dynamics.Internal.Dynamics
        (-- * Dynamics
         Dynamics(..),
+        DynamicsLift(..),
         Point(..),
-        Specs(..),
-        Method(..),
-        Run(..),
-        runDynamics1,
-        runDynamics1_,
+        runDynamicsInStart,
+        runDynamicsInFinal,
         runDynamics,
-        runDynamics_,
-        runDynamicsIO,
-        runDynamicsSeries1,
-        runDynamicsSeries1_,
-        runDynamicsSeries,
-        runDynamicsSeries_,
-        printDynamics1,
-        printDynamics,
         -- * Utilities
         basicTime,
         iterationBnds,
@@ -46,6 +30,8 @@
 import Control.Monad
 import Control.Monad.Trans
 
+import Simulation.Aivika.Dynamics.Internal.Simulation
+
 --
 -- The Dynamics Monad
 --
@@ -65,24 +51,6 @@
                      pointIteration :: Int,  -- ^ the current iteration
                      pointPhase :: Int       -- ^ the current phase
                    } deriving (Eq, Ord, Show)
-
--- | It defines the simulation specs.
-data Specs = Specs { spcStartTime :: Double,    -- ^ the start time
-                     spcStopTime :: Double,     -- ^ the stop time
-                     spcDT :: Double,           -- ^ the integration time step
-                     spcMethod :: Method        -- ^ the integration method
-                   } deriving (Eq, Ord, Show)
-
--- | It defines the integration method.
-data Method = Euler          -- ^ Euler's method
-            | RungeKutta2    -- ^ the 2nd order Runge-Kutta method
-            | RungeKutta4    -- ^ the 4th order Runge-Kutta method
-            deriving (Eq, Ord, Show)
-
--- | It defined the simulation run as part of some experiment.
-data Run = Run { runIndex :: Int,    -- ^ the current simulation run
-                 runCount :: Int     -- ^ the total number of runs in this experiment
-               } deriving (Eq, Ord, Show)
            
 -- | Returns the iterations starting from zero.
 iterations :: Specs -> [Int]
@@ -164,174 +132,44 @@
      let Dynamics m' = k a
      m' p
 
-subrunDynamics1 :: Dynamics a -> Specs -> Run -> IO a
-subrunDynamics1 (Dynamics m) sc r =
-  do let n = iterationHiBnd sc
-         t = basicTime sc n 0
+-- | Run the dynamic process in the initial simulation point.
+runDynamicsInStart :: Dynamics a -> Simulation a
+runDynamicsInStart (Dynamics m) =
+  Simulation $ \r ->
+  do let sc = runSpecs r 
+         n  = 0
+         t  = spcStartTime sc
      m Point { pointSpecs = sc,
                pointRun = r,
                pointTime = t,
                pointIteration = n,
                pointPhase = 0 }
 
-subrunDynamics1_ :: Dynamics a -> Specs -> Run -> IO ()
-subrunDynamics1_ (Dynamics m) sc r =
-  do let n = iterationHiBnd sc
-         t = basicTime sc n 0
+-- | Run the dynamic process in the final simulation point.
+runDynamicsInFinal :: Dynamics a -> Simulation a
+runDynamicsInFinal (Dynamics m) =
+  Simulation $ \r ->
+  do let sc = runSpecs r 
+         n  = iterationHiBnd sc
+         t  = basicTime sc n 0
      m Point { pointSpecs = sc,
                pointRun = r,
                pointTime = t,
                pointIteration = n,
                pointPhase = 0 }
-     return ()
 
-subrunDynamics :: Dynamics a -> Specs -> Run -> [IO a]
-subrunDynamics (Dynamics m) sc r =
-  do let (nl, nu) = iterationBnds sc
-         point n = Point { pointSpecs = sc,
-                           pointRun = r,
-                           pointTime = basicTime sc n 0,
-                           pointIteration = n,
-                           pointPhase = 0 }
-     map (m . point) [nl .. nu]
-
-subrunDynamics_ :: Dynamics a -> Specs -> Run -> IO ()
-subrunDynamics_ (Dynamics m) sc r =
-  do let (nl, nu) = iterationBnds sc
+-- | Run the dynamic process in all integration time points
+runDynamics :: Dynamics a -> Simulation [IO a]
+runDynamics (Dynamics m) =
+  Simulation $ \r ->
+  do let sc = runSpecs r
+         (nl, nu) = iterationBnds sc
          point n = Point { pointSpecs = sc,
                            pointRun = r,
                            pointTime = basicTime sc n 0,
                            pointIteration = n,
                            pointPhase = 0 }
-     mapM_ (m . point) [nl .. nu]
-
--- | Run the simulation and return the result in the last 
--- time point using the specified simulation specs.
-runDynamics1 :: Dynamics (Dynamics a) -> Specs -> IO a
-runDynamics1 (Dynamics m) sc = 
-  do let r = Run { runIndex = 1, runCount = 1 }
-     d <- m Point { pointSpecs = sc,
-                    pointRun = r,
-                    pointTime = spcStartTime sc,
-                    pointIteration = 0,
-                    pointPhase = 0 }
-     subrunDynamics1 d sc r
-
--- | Run the simulation and return the result in the last 
--- time point using the specified simulation specs.
-runDynamics1_ :: Dynamics (Dynamics a) -> Specs -> IO ()
-runDynamics1_ (Dynamics m) sc = 
-  do let r = Run { runIndex = 1, runCount = 1 }
-     d <- m Point { pointSpecs = sc,
-                    pointRun = r,
-                    pointTime = spcStartTime sc,
-                    pointIteration = 0,
-                    pointPhase = 0 }
-     subrunDynamics1_ d sc r
-
--- | Run the simulation and return the results in all 
--- integration time points using the specified simulation specs.
-runDynamics :: Dynamics (Dynamics a) -> Specs -> IO [a]
-runDynamics (Dynamics m) sc = 
-  do let r = Run { runIndex = 1, runCount = 1 }
-     d <- m Point { pointSpecs = sc,
-                    pointRun = r,
-                    pointTime = spcStartTime sc,
-                    pointIteration = 0,
-                    pointPhase = 0 }
-     sequence $ subrunDynamics d sc r
-
--- | Run the simulation and return the results in all 
--- integration time points using the specified simulation specs.
-runDynamics_ :: Dynamics (Dynamics a) -> Specs -> IO ()
-runDynamics_ (Dynamics m) sc = 
-  do let r = Run { runIndex = 1, runCount = 1 }
-     d <- m Point { pointSpecs = sc,
-                    pointRun = r,
-                    pointTime = spcStartTime sc,
-                    pointIteration = 0,
-                    pointPhase = 0 }
-     sequence_ $ subrunDynamics d sc r
-
--- | Run the simulation and return the results in all 
--- integration time points using the specified simulation specs.
-runDynamicsIO :: Dynamics (Dynamics a) -> Specs -> IO [IO a]
-runDynamicsIO (Dynamics m) sc =
-  do let r = Run { runIndex = 1, runCount = 1 }
-     d <- m Point { pointSpecs = sc,
-                    pointRun = r,
-                    pointTime = spcStartTime sc,
-                    pointIteration = 0,
-                    pointPhase = 0 }
-     return $ subrunDynamics d sc r
-
--- | Run an experiment consisting of the given number of simulations, where each 
--- model is created and then requested in the last integration time point using 
--- the specified specs.
-runDynamicsSeries1_ :: Dynamics (Dynamics a) -> Specs -> Int -> [IO ()]
-runDynamicsSeries1_ (Dynamics m) sc runs = map f [1 .. runs]
-  where f i =
-          do let r = Run { runIndex = i, runCount = runs }
-             d <- m Point { pointSpecs = sc,
-                            pointRun = r,
-                            pointTime = spcStartTime sc,
-                            pointIteration = 0,
-                            pointPhase = 0 }
-             subrunDynamics1_ d sc r
-
--- | Run an experiment consisting of the given number of simulations, where each 
--- model is created and then requested sequentially in all integration time points 
--- using the specified specs.
-runDynamicsSeries_ :: Dynamics (Dynamics a) -> Specs -> Int -> [IO ()]
-runDynamicsSeries_ (Dynamics m) sc runs = map f [1 .. runs]
-  where f i =
-          do let r = Run { runIndex = i, runCount = runs }
-             d <- m Point { pointSpecs = sc,
-                            pointRun = r,
-                            pointTime = spcStartTime sc,
-                            pointIteration = 0,
-                            pointPhase = 0 }
-             subrunDynamics_ d sc r
-
--- | Run an experiment consisting of the given number of simulations, where each 
--- model is created and then requested in the last integration time point using 
--- the specified specs.
-runDynamicsSeries1 :: Dynamics (Dynamics a) -> Specs -> Int -> [IO a]
-runDynamicsSeries1 (Dynamics m) sc runs = map f [1 .. runs]
-  where f i =
-          do let r = Run { runIndex = i, runCount = runs }
-             d <- m Point { pointSpecs = sc,
-                            pointRun = r,
-                            pointTime = spcStartTime sc,
-                            pointIteration = 0,
-                            pointPhase = 0 }
-             subrunDynamics1 d sc r
-
--- | Run an experiment consisting of the given number of simulations, where each 
--- model is created and then requested sequentially in all integration time points 
--- using the specified specs.
-runDynamicsSeries :: Dynamics (Dynamics a) -> Specs -> Int -> [IO [a]]
-runDynamicsSeries (Dynamics m) sc runs = map f [1 .. runs]
-  where f i =
-          do let r = Run { runIndex = i, runCount = runs }
-             d <- m Point { pointSpecs = sc,
-                            pointRun = r,
-                            pointTime = spcStartTime sc,
-                            pointIteration = 0,
-                            pointPhase = 0 }
-             sequence $ subrunDynamics d sc r
-
--- | Run the simulation and print the result in the last 
--- time point using the specified simulation specs.
-printDynamics1 :: (Show a) => Dynamics (Dynamics a) -> Specs -> IO ()
-printDynamics1 m sc = runDynamics1 m sc >>= print
-
--- | Run the simulation and print lazily the results in all
--- integration time points using the specified simulation specs.
-printDynamics :: (Show a) => Dynamics (Dynamics a) -> Specs -> IO ()
-printDynamics m sc = runDynamicsIO m sc >>= loop
-  where loop [] = return ()
-        loop (x : xs) = do { a <- x; print a; loop xs }
+     return $ map (m . point) [nl .. nu]
 
 instance Functor Dynamics where
   fmap = liftMD
@@ -387,3 +225,17 @@
 
 instance MonadIO Dynamics where
   liftIO m = Dynamics $ const m
+
+instance SimulationLift Dynamics where
+  liftSimulation = liftDS
+    
+liftDS :: Simulation a -> Dynamics a
+{-# INLINE liftDS #-}
+liftDS (Simulation m) =
+  Dynamics $ \p -> m $ pointRun p
+
+-- | A type class to lift the 'Dynamics' computations in other monads.
+class Monad m => DynamicsLift m where
+  
+  -- | Lift the specified 'Dynamics' computation in another monad.
+  liftDynamics :: Dynamics a -> m a
diff --git a/Simulation/Aivika/Dynamics/Internal/Fold.hs b/Simulation/Aivika/Dynamics/Internal/Fold.hs
--- a/Simulation/Aivika/Dynamics/Internal/Fold.hs
+++ b/Simulation/Aivika/Dynamics/Internal/Fold.hs
@@ -11,14 +11,15 @@
 -- any dynamic process in the integration time points.
 --
 module Simulation.Aivika.Dynamics.Internal.Fold
-       (foldD1,
-        foldD,
-        divideD) where
+       (foldDynamics1,
+        foldDynamics,
+        divideDynamics) where
 
 import Data.IORef
 import Control.Monad
 import Control.Monad.Trans
 
+import Simulation.Aivika.Dynamics.Internal.Simulation
 import Simulation.Aivika.Dynamics.Internal.Dynamics
 import Simulation.Aivika.Dynamics.Internal.Interpolate
 import Simulation.Aivika.Dynamics.Internal.Memo
@@ -31,11 +32,11 @@
 -- the integration time points. The accumulator values are transformed
 -- according to the first argument, which should be either function 
 -- 'memo0' or 'umemo0'.
-foldD1 :: (Dynamics a -> Dynamics (Dynamics a))
-         -> (a -> a -> a) 
-         -> Dynamics a 
-         -> Dynamics (Dynamics a)
-foldD1 tr f (Dynamics m) =
+foldDynamics1 :: (Dynamics a -> Simulation (Dynamics a))
+                 -> (a -> a -> a) 
+                 -> Dynamics a 
+                 -> Simulation (Dynamics a)
+foldDynamics1 tr f (Dynamics m) =
   do r <- liftIO $ newIORef m
      let z = Dynamics $ \p ->
            case pointIteration p of
@@ -57,12 +58,12 @@
 -- the integration time points. The accumulator values are transformed
 -- according to the first argument, which should be either function
 -- 'memo0' or 'umemo0'.
-foldD :: (Dynamics a -> Dynamics (Dynamics a))
-        -> (a -> b -> a) 
-        -> a
-        -> Dynamics b 
-        -> Dynamics (Dynamics a)
-foldD tr f acc (Dynamics m) =
+foldDynamics :: (Dynamics a -> Simulation (Dynamics a))
+                -> (a -> b -> a) 
+                -> a
+                -> Dynamics b 
+                -> Simulation (Dynamics a)
+foldDynamics tr f acc (Dynamics m) =
   do r <- liftIO $ newIORef $ const $ return acc
      let z = Dynamics $ \p ->
            case pointIteration p of
@@ -84,8 +85,8 @@
 -- | Divide the values in integration time points by the number of
 -- the current iteration. It can be useful for statistic functions in
 -- combination with the fold.
-divideD :: Dynamics Double -> Dynamics Double
-divideD (Dynamics m) = 
+divideDynamics :: Dynamics Double -> Dynamics Double
+divideDynamics (Dynamics m) = 
   discrete $ Dynamics $ \p ->
   do a <- m p
      return $ a / fromInteger (toInteger (pointIteration p + 1))
diff --git a/Simulation/Aivika/Dynamics/Internal/Interpolate.hs b/Simulation/Aivika/Dynamics/Internal/Interpolate.hs
--- a/Simulation/Aivika/Dynamics/Internal/Interpolate.hs
+++ b/Simulation/Aivika/Dynamics/Internal/Interpolate.hs
@@ -14,16 +14,16 @@
 --
 
 module Simulation.Aivika.Dynamics.Internal.Interpolate
-       (initD,
+       (initDynamics,
         discrete,
         interpolate) where
 
 import Simulation.Aivika.Dynamics.Internal.Dynamics
 
 -- | Return the initial value.
-initD :: Dynamics a -> Dynamics a
-{-# INLINE initD #-}
-initD (Dynamics m) =
+initDynamics :: Dynamics a -> Dynamics a
+{-# INLINE initDynamics #-}
+initDynamics (Dynamics m) =
   Dynamics $ \p ->
   if pointIteration p == 0 && pointPhase p == 0 then
     m p
diff --git a/Simulation/Aivika/Dynamics/Internal/Memo.hs b/Simulation/Aivika/Dynamics/Internal/Memo.hs
--- a/Simulation/Aivika/Dynamics/Internal/Memo.hs
+++ b/Simulation/Aivika/Dynamics/Internal/Memo.hs
@@ -19,13 +19,14 @@
         umemo,
         memo0,
         umemo0,
-        iterateD) where
+        iterateDynamics) where
 
 import Data.Array
 import Data.Array.IO
 import Data.IORef
 import Control.Monad
 
+import Simulation.Aivika.Dynamics.Internal.Simulation
 import Simulation.Aivika.Dynamics.Internal.Dynamics
 import Simulation.Aivika.Dynamics.Internal.Interpolate
 
@@ -37,11 +38,11 @@
 
 -- | Memoize and order the computation in the integration time points using 
 -- the interpolation that knows of the Runge-Kutta method.
-memo :: Dynamics e -> Dynamics (Dynamics e)
+memo :: Dynamics e -> Simulation (Dynamics e)
 {-# INLINE memo #-}
 memo (Dynamics m) = 
-  Dynamics $ \p ->
-  do let sc = pointSpecs p
+  Simulation $ \r ->
+  do let sc = runSpecs r
          (phl, phu) = phaseBnds sc
          (nl, nu)   = iterationBnds sc
      arr   <- newMemoArray_ ((phl, nl), (phu, nu))
@@ -74,11 +75,11 @@
 
 -- | This is a more efficient version the 'memo' function which uses 
 -- an unboxed array to store the values.
-umemo :: (MArray IOUArray e IO) => Dynamics e -> Dynamics (Dynamics e)
+umemo :: (MArray IOUArray e IO) => Dynamics e -> Simulation (Dynamics e)
 {-# INLINE umemo #-}
 umemo (Dynamics m) = 
-  Dynamics $ \p ->
-  do let sc = pointSpecs p
+  Simulation $ \r ->
+  do let sc = runSpecs r
          (phl, phu) = phaseBnds sc
          (nl, nu)   = iterationBnds sc
      arr   <- newMemoUArray_ ((phl, nl), (phu, nu))
@@ -116,11 +117,11 @@
 -- difference when we request for values in the intermediate time points
 -- that are used by this method to integrate. In general case you should 
 -- prefer the 'memo0' function above 'memo'.
-memo0 :: Dynamics e -> Dynamics (Dynamics e)
+memo0 :: Dynamics e -> Simulation (Dynamics e)
 {-# INLINE memo0 #-}
 memo0 (Dynamics m) = 
-  Dynamics $ \p ->
-  do let sc   = pointSpecs p
+  Simulation $ \r ->
+  do let sc   = runSpecs r
          bnds = iterationBnds sc
      arr  <- newMemoArray_ bnds
      nref <- newIORef 0
@@ -144,11 +145,11 @@
 
 -- | This is a more efficient version the 'memo0' function which uses 
 -- an unboxed array to store the values.
-umemo0 :: (MArray IOUArray e IO) => Dynamics e -> Dynamics (Dynamics e)
+umemo0 :: (MArray IOUArray e IO) => Dynamics e -> Simulation (Dynamics e)
 {-# INLINE umemo0 #-}
 umemo0 (Dynamics m) = 
-  Dynamics $ \p ->
-  do let sc   = pointSpecs p
+  Simulation $ \r ->
+  do let sc   = runSpecs r
          bnds = iterationBnds sc
      arr  <- newMemoUArray_ bnds
      nref <- newIORef 0
@@ -174,11 +175,11 @@
 -- the integration time points. It is equivalent to a call of the
 -- 'memo0' function but significantly more efficient, for the array 
 -- is not created.
-iterateD :: Dynamics () -> Dynamics (Dynamics ())
-{-# INLINE iterateD #-}
-iterateD (Dynamics m) = 
-  Dynamics $ \p ->
-  do let sc = pointSpecs p
+iterateDynamics :: Dynamics () -> Simulation (Dynamics ())
+{-# INLINE iterateDynamics #-}
+iterateDynamics (Dynamics m) = 
+  Simulation $ \r ->
+  do let sc = runSpecs r
      nref <- newIORef 0
      let r p =
            do let sc = pointSpecs p
diff --git a/Simulation/Aivika/Dynamics/Internal/Process.hs b/Simulation/Aivika/Dynamics/Internal/Process.hs
--- a/Simulation/Aivika/Dynamics/Internal/Process.hs
+++ b/Simulation/Aivika/Dynamics/Internal/Process.hs
@@ -1,16 +1,19 @@
 
 -- |
 -- Module     : Simulation.Aivika.Dynamics.Internal.Process
--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>
 -- License    : BSD3
 -- Maintainer : David Sorokin <david.sorokin@gmail.com>
 -- Stability  : experimental
 -- Tested with: GHC 7.0.3
 --
 -- A value in the 'Process' monad represents a discontinuous process that 
--- can suspend and resume at any time. It behaves like a dynamic process too. 
--- Any value in the 'Dynamics' monad can be lifted to the Process monad. 
--- Moreover, a value in the Process monad can be run in the Dynamics monad.
+-- can suspend in any simulation time point and then resume later in the same 
+-- or another time point. 
+-- 
+-- The process of this type behaves like a dynamic process too. So, any value 
+-- in the 'Dynamics' monad can be lifted to the Process monad. Moreover, 
+-- a value in the Process monad can be run in the Dynamics monad.
 --
 -- A value of the 'ProcessID' type is just an identifier of such a process.
 --
@@ -26,20 +29,21 @@
         processID,
         runProcess) where
 
+import Data.Maybe
 import Data.IORef
 import Control.Monad
 import Control.Monad.Trans
 
+import Simulation.Aivika.Dynamics.Internal.Simulation
 import Simulation.Aivika.Dynamics.Internal.Dynamics
 import Simulation.Aivika.Dynamics.Internal.Cont
-import Simulation.Aivika.Dynamics.Lift
 import Simulation.Aivika.Dynamics.EventQueue
 
 -- | Represents a process identificator.
 data ProcessID = 
   ProcessID { processQueue   :: EventQueue,  -- ^ Return the event queue.
               processStarted :: IORef Bool,
-              processCont    :: IORef (Maybe (Dynamics (() -> IO ()))) }
+              processCont    :: IORef (Maybe (() -> Dynamics ())) }
 
 -- | Specifies a discontinuous process that can suspend at any time
 -- and then resume later.
@@ -50,9 +54,9 @@
 holdProcess dt =
   Process $ \pid ->
   Cont $ \c ->
-  Dynamics $ \ps ->
-  do let Dynamics m = enqueueCont (processQueue pid) (pointTime ps + dt) c
-     m ps
+  Dynamics $ \p ->
+  do let Dynamics m = enqueueCont (processQueue pid) (pointTime p + dt) c
+     m p
 
 -- | Passivate the process.
 passivateProcess :: Process ()
@@ -67,41 +71,35 @@
        Just _  -> error "Cannot passivate the process twice: passivate"
 
 -- | Test whether the process with the specified ID is passivated.
-processPassive :: ProcessID -> Process Bool
+processPassive :: ProcessID -> Dynamics Bool
 processPassive pid =
-  Process $ \_ ->
-  Cont $ \(Dynamics c) ->
   Dynamics $ \p ->
-  do cont' <- c p
+  do let Dynamics m = queueRun $ processQueue pid
+     m p
      let x = processCont pid
      a <- readIORef x
-     case a of
-       Nothing -> cont' False
-       Just _  -> cont' True
+     return $ isJust a
 
 -- | Reactivate a process with the specified ID.
-reactivateProcess :: ProcessID -> Process ()
+reactivateProcess :: ProcessID -> Dynamics ()
 reactivateProcess pid =
-  Process $ \pid' ->
-  Cont $ \c@(Dynamics cont) ->
   Dynamics $ \p ->
-  do let x = processCont pid
+  do let Dynamics m = queueRun $ processQueue pid
+     m p
+     let x = processCont pid
      a <- readIORef x
      case a of
-       Nothing ->
-         do cont' <- cont p
-            cont' ()
-       Just (Dynamics cont2) ->
+       Nothing -> 
+         return ()
+       Just c ->
          do writeIORef x Nothing
-            let Dynamics m = enqueueCont (processQueue pid') (pointTime p) c
+            let Dynamics m  = enqueueCont (processQueue pid) (pointTime p) c
             m p
-            cont2' <- cont2 p
-            cont2' ()
 
 -- | Start the process with the specified ID at the desired time.
 runProcess :: Process () -> ProcessID -> Double -> Dynamics ()
 runProcess (Process p) pid t =
-  runCont m r
+  runCont m return
     where m = do y <- liftIO $ readIORef (processStarted pid)
                  if y 
                    then error $
@@ -110,14 +108,13 @@
                    else liftIO $ writeIORef (processStarted pid) True
                  Cont $ \c -> enqueueCont (processQueue pid) t c
                  p pid
-          r = let f () = return () in return f
 
 -- | Return the current process ID.
 processID :: Process ProcessID
 processID = Process $ \pid -> return pid
 
 -- | Create a new process ID.
-newProcessID :: EventQueue -> Dynamics ProcessID
+newProcessID :: EventQueue -> Simulation ProcessID
 newProcessID q =
   do x <- liftIO $ newIORef Nothing
      y <- liftIO $ newIORef False
@@ -135,9 +132,12 @@
 instance Functor Process where
   fmap = liftM
 
-instance Lift Process where
-  liftD = liftP
+instance SimulationLift Process where
+  liftSimulation = liftSP
   
+instance DynamicsLift Process where
+  liftDynamics = liftDP
+  
 instance MonadIO Process where
   liftIO = liftIOP
   
@@ -153,9 +153,13 @@
      let Process m' = k a
      m' pid
 
-liftP :: Dynamics a -> Process a
-{-# INLINE liftP #-}
-liftP m = Process $ \pid -> liftD m
+liftSP :: Simulation a -> Process a
+{-# INLINE liftSP #-}
+liftSP m = Process $ \pid -> liftSimulation m
+
+liftDP :: Dynamics a -> Process a
+{-# INLINE liftDP #-}
+liftDP m = Process $ \pid -> liftDynamics m
 
 liftIOP :: IO a -> Process a
 {-# INLINE liftIOP #-}
diff --git a/Simulation/Aivika/Dynamics/Internal/Simulation.hs b/Simulation/Aivika/Dynamics/Internal/Simulation.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Internal/Simulation.hs
@@ -0,0 +1,113 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Internal.Simulation
+-- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- The module defines the 'Simulation' monad that represents a simulation run.
+-- 
+module Simulation.Aivika.Dynamics.Internal.Simulation
+       (-- * Simulation
+        Simulation(..),
+        SimulationLift(..),
+        Specs(..),
+        Method(..),
+        Run(..),
+        runSimulation,
+        runSimulations,
+        -- * Utilities
+        simulationIndex,
+        simulationCount,
+        simulationSpecs) where
+
+import Control.Monad
+import Control.Monad.Trans
+
+--
+-- The Simulation Monad
+--
+-- A value of the Simulation monad represents a simulation run.
+--
+
+-- | A value in the 'Simulation' monad represents a simulation run.
+newtype Simulation a = Simulation (Run -> IO a)
+
+-- | It defines the simulation specs.
+data Specs = Specs { spcStartTime :: Double,    -- ^ the start time
+                     spcStopTime :: Double,     -- ^ the stop time
+                     spcDT :: Double,           -- ^ the integration time step
+                     spcMethod :: Method        -- ^ the integration method
+                   } deriving (Eq, Ord, Show)
+
+-- | It defines the integration method.
+data Method = Euler          -- ^ Euler's method
+            | RungeKutta2    -- ^ the 2nd order Runge-Kutta method
+            | RungeKutta4    -- ^ the 4th order Runge-Kutta method
+            deriving (Eq, Ord, Show)
+
+-- | It indentifies the simulation run.
+data Run = Run { runSpecs :: Specs,  -- ^ the simulation specs
+                 runIndex :: Int,    -- ^ the current simulation run index
+                 runCount :: Int     -- ^ the total number of runs in this experiment
+               } deriving (Eq, Ord, Show)
+
+instance Monad Simulation where
+  return  = returnS
+  m >>= k = bindS m k
+
+returnS :: a -> Simulation a
+returnS a = Simulation (\r -> return a)
+
+bindS :: Simulation a -> (a -> Simulation b) -> Simulation b
+bindS (Simulation m) k = 
+  Simulation $ \r -> 
+  do a <- m r
+     let Simulation m' = k a
+     m' r
+
+-- | Run the simulation using the specified specs.
+runSimulation :: Simulation a -> Specs -> IO a
+runSimulation (Simulation m) sc =
+  m Run { runSpecs = sc,
+          runIndex = 1,
+          runCount = 1 }
+
+-- | Run the given number of simulations using the specified specs, 
+--   where each simulation is distinguished by its index 'simulationIndex'.
+runSimulations :: Simulation a -> Specs -> Int -> [IO a]
+runSimulations (Simulation m) sc runs = map f [1 .. runs]
+  where f i = m Run { runSpecs = sc,
+                      runIndex = i,
+                      runCount = runs }
+
+-- | Return the run index for the current simulation.
+simulationIndex :: Simulation Int
+simulationIndex = Simulation $ return . runIndex
+
+-- | Return the number of simulations currently run.
+simulationCount :: Simulation Int
+simulationCount = Simulation $ return . runCount
+
+-- | Return the simulation specs
+simulationSpecs :: Simulation Specs
+simulationSpecs = Simulation $ return . runSpecs
+
+instance Functor Simulation where
+  fmap = liftMS
+
+liftMS :: (a -> b) -> Simulation a -> Simulation b
+{-# INLINE liftMS #-}
+liftMS f (Simulation x) =
+  Simulation $ \r -> do { a <- x r; return $ f a }
+
+instance MonadIO Simulation where
+  liftIO m = Simulation $ const m
+
+-- | A type class to lift the simulation computations in other monads.
+class Monad m => SimulationLift m where
+  
+  -- | Lift the specified 'Simulation' computation in another monad.
+  liftSimulation :: Simulation a -> m a
diff --git a/Simulation/Aivika/Dynamics/Internal/Time.hs b/Simulation/Aivika/Dynamics/Internal/Time.hs
--- a/Simulation/Aivika/Dynamics/Internal/Time.hs
+++ b/Simulation/Aivika/Dynamics/Internal/Time.hs
@@ -11,11 +11,12 @@
 --
 
 module Simulation.Aivika.Dynamics.Internal.Time
-       (starttime,
-        stoptime,
-        dt,
+       (starttime, 
+        stoptime, 
+        dt, 
         time) where
 
+import Simulation.Aivika.Dynamics.Internal.Simulation
 import Simulation.Aivika.Dynamics.Internal.Dynamics
 
 -- | Return the start simulation time.
diff --git a/Simulation/Aivika/Dynamics/Lift.hs b/Simulation/Aivika/Dynamics/Lift.hs
deleted file mode 100644
--- a/Simulation/Aivika/Dynamics/Lift.hs
+++ /dev/null
@@ -1,21 +0,0 @@
-
--- |
--- Module     : Simulation.Aivika.Dynamics.Lift
--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
--- License    : BSD3
--- Maintainer : David Sorokin <david.sorokin@gmail.com>
--- Stability  : experimental
--- Tested with: GHC 7.0.3
---
--- This module defines the 'liftD' function that allows embedding
--- the 'Dynamics' computation.
---
-module Simulation.Aivika.Dynamics.Lift (Lift(..)) where
-
-import Simulation.Aivika.Dynamics
-
--- | The 'Lift' class defines a type which the 'Dynamics' 
--- computation can be lifted to.
-class Lift m where
-  -- | Lift the computation.
-  liftD :: Dynamics a -> m a
diff --git a/Simulation/Aivika/Dynamics/Parameter.hs b/Simulation/Aivika/Dynamics/Parameter.hs
--- a/Simulation/Aivika/Dynamics/Parameter.hs
+++ b/Simulation/Aivika/Dynamics/Parameter.hs
@@ -20,6 +20,7 @@
 import qualified Data.Map as M
 import Control.Concurrent.MVar
 
+import Simulation.Aivika.Dynamics.Internal.Simulation
 import Simulation.Aivika.Dynamics.Internal.Dynamics
 
 -- | Create a thread-safe parameter that returns always the same value during the simulation run, 
diff --git a/Simulation/Aivika/Dynamics/Process.hs b/Simulation/Aivika/Dynamics/Process.hs
--- a/Simulation/Aivika/Dynamics/Process.hs
+++ b/Simulation/Aivika/Dynamics/Process.hs
@@ -8,9 +8,12 @@
 -- Tested with: GHC 7.0.3
 --
 -- A value in the 'Process' monad represents a discontinuous process that 
--- can suspend and resume at any time. It behaves like a dynamic process too. 
--- Any value in the 'Dynamics' monad can be lifted to the Process monad. 
--- Moreover, a value in the Process monad can be run in the Dynamics monad.
+-- can suspend in any simulation time point and then resume later in the same 
+-- or another time point. 
+-- 
+-- The process of this type behaves like a dynamic process too. So, any value 
+-- in the 'Dynamics' monad can be lifted to the Process monad. Moreover, 
+-- a value in the Process monad can be run in the Dynamics monad.
 --
 -- A value of the 'ProcessID' type is just an identifier of such a process.
 --
diff --git a/Simulation/Aivika/Dynamics/Random.hs b/Simulation/Aivika/Dynamics/Random.hs
--- a/Simulation/Aivika/Dynamics/Random.hs
+++ b/Simulation/Aivika/Dynamics/Random.hs
@@ -15,22 +15,23 @@
 module Simulation.Aivika.Dynamics.Random 
        (newRandom, newNormal, normalGen) where
 
-import Random
+import System.Random
 import Data.IORef
 import Control.Monad.Trans
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Simulation
 import Simulation.Aivika.Dynamics.Base
 
 -- | Return the uniform random numbers between 0.0 and 1.0 in
 -- the integration time points.
-newRandom :: Dynamics (Dynamics Double)
+newRandom :: Simulation (Dynamics Double)
 newRandom =
   memo0 $ liftIO $ getStdRandom random
      
 -- | Return the normal random numbers with mean 0.0 and variance 1.0 in
 -- the integration time points.
-newNormal :: Dynamics (Dynamics Double)
+newNormal :: Simulation (Dynamics Double)
 newNormal =
   do g <- liftIO normalGen
      memo0 $ liftIO g
diff --git a/Simulation/Aivika/Dynamics/Ref.hs b/Simulation/Aivika/Dynamics/Ref.hs
--- a/Simulation/Aivika/Dynamics/Ref.hs
+++ b/Simulation/Aivika/Dynamics/Ref.hs
@@ -11,15 +11,17 @@
 --
 module Simulation.Aivika.Dynamics.Ref
        (Ref,
-        newRef,
         refQueue,
+        newRef,
         readRef,
         writeRef,
         modifyRef) where
 
 import Data.IORef
+import Control.Monad
 import Control.Monad.Trans
 
+import Simulation.Aivika.Dynamics.Internal.Simulation
 import Simulation.Aivika.Dynamics.Internal.Dynamics
 import Simulation.Aivika.Dynamics.EventQueue
 
@@ -32,7 +34,7 @@
         refValue :: IORef a }
 
 -- | Create a new reference bound to the specified event queue.
-newRef :: EventQueue -> a -> Dynamics (Ref a)
+newRef :: EventQueue -> a -> Simulation (Ref a)
 newRef q a =
   do x <- liftIO $ newIORef a
      return Ref { refQueue = q,
diff --git a/Simulation/Aivika/Dynamics/Resource.hs b/Simulation/Aivika/Dynamics/Resource.hs
--- a/Simulation/Aivika/Dynamics/Resource.hs
+++ b/Simulation/Aivika/Dynamics/Resource.hs
@@ -1,14 +1,14 @@
 
 -- |
 -- Module     : Simulation.Aivika.Dynamics.Resource
--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>
 -- License    : BSD3
 -- Maintainer : David Sorokin <david.sorokin@gmail.com>
 -- Stability  : experimental
 -- Tested with: GHC 7.0.3
 --
 -- This module defines a limited resource which can be acquired and 
--- then released by the discontinuous process 'DynamicProc'.
+-- then released by the discontinuous process 'Process'.
 --
 module Simulation.Aivika.Dynamics.Resource
        (Resource,
@@ -22,6 +22,7 @@
 import Data.IORef
 import Control.Monad
 
+import Simulation.Aivika.Dynamics.Internal.Simulation
 import Simulation.Aivika.Dynamics.Internal.Dynamics
 import Simulation.Aivika.Dynamics.Internal.Cont
 import Simulation.Aivika.Dynamics.Internal.Process
@@ -35,15 +36,15 @@
              resourceInitCount :: Int,
              -- ^ Return the initial count of the resource.
              resourceCountRef  :: IORef Int, 
-             resourceWaitQueue :: Q.Queue (Dynamics (() -> IO ()))}
+             resourceWaitQueue :: Q.Queue (() -> Dynamics ())}
 
 instance Eq Resource where
   x == y = resourceCountRef x == resourceCountRef y  -- unique references
 
 -- | Create a new resource with the specified initial count.
-newResource :: EventQueue -> Int -> Dynamics Resource
+newResource :: EventQueue -> Int -> Simulation Resource
 newResource q initCount =
-  Dynamics $ \p ->
+  Simulation $ \r ->
   do countRef  <- newIORef initCount
      waitQueue <- Q.newQueue
      return Resource { resourceQueue     = q,
@@ -52,14 +53,12 @@
                        resourceWaitQueue = waitQueue }
 
 -- | Return the current count of the resource.
-resourceCount :: Resource -> Process Int
+resourceCount :: Resource -> Dynamics Int
 resourceCount r =
-  Process $ \_ ->
-  Cont $ \(Dynamics c) ->
   Dynamics $ \p ->
-  do cont' <- c p 
-     a <- readIORef (resourceCountRef r)
-     cont' a
+  do let Dynamics m = queueRun (resourceQueue r)
+     m p
+     readIORef (resourceCountRef r)
 
 -- | Request for the resource decreasing its count in case of success,
 -- otherwise suspending the discontinuous process until some other 
@@ -67,22 +66,22 @@
 requestResource :: Resource -> Process ()
 requestResource r =
   Process $ \_ ->
-  Cont $ \c@(Dynamics cont) ->
+  Cont $ \c ->
   Dynamics $ \p ->
   do a <- readIORef (resourceCountRef r)
      if a == 0 
        then Q.enqueue (resourceWaitQueue r) c
        else do let a' = a - 1
                a' `seq` writeIORef (resourceCountRef r) a'
-               cont' <- cont p
-               cont' ()
+               let Dynamics m = c ()
+               m p
 
 -- | Release the resource increasing its count and resuming one of the
 -- previously suspended processes as possible.
 releaseResource :: Resource -> Process ()
 releaseResource r =
   Process $ \_ ->
-  Cont $ \(Dynamics c) ->
+  Cont $ \c ->
   Dynamics $ \p ->
   do a <- readIORef (resourceCountRef r)
      let a' = a + 1
@@ -97,5 +96,5 @@
                Q.dequeue (resourceWaitQueue r)
                let Dynamics m = enqueueCont (resourceQueue r) (pointTime p) c2
                m p
-     cont' <- c p
-     cont' ()
+     let Dynamics m' = c ()
+     m' p
diff --git a/Simulation/Aivika/Dynamics/Simulation.hs b/Simulation/Aivika/Dynamics/Simulation.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Simulation.hs
@@ -0,0 +1,23 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Simulation
+-- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- The module defines the 'Simulation' monad representing a simulation run.
+--
+module Simulation.Aivika.Dynamics.Simulation
+       (Simulation,
+        SimulationLift(..),
+        Specs(..),
+        Method(..),
+        runSimulation,
+        runSimulations,
+        simulationIndex,
+        simulationCount,
+        simulationSpecs) where
+
+import Simulation.Aivika.Dynamics.Internal.Simulation
diff --git a/Simulation/Aivika/Dynamics/SystemDynamics.hs b/Simulation/Aivika/Dynamics/SystemDynamics.hs
--- a/Simulation/Aivika/Dynamics/SystemDynamics.hs
+++ b/Simulation/Aivika/Dynamics/SystemDynamics.hs
@@ -1,5 +1,5 @@
 
-{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleContexts, BangPatterns #-}
 
 -- |
 -- Module     : Simulation.Aivika.Dynamics.SystemDynamics
@@ -14,8 +14,8 @@
 
 module Simulation.Aivika.Dynamics.SystemDynamics
        (-- * Maximum and Minimum
-        maxD,
-        minD,
+        maxDynamics,
+        minDynamics,
         -- * Integrals
         Integ,
         newInteg,
@@ -40,6 +40,7 @@
 import Control.Monad
 import Control.Monad.Trans
 
+import Simulation.Aivika.Dynamics.Internal.Simulation
 import Simulation.Aivika.Dynamics.Internal.Dynamics
 import Simulation.Aivika.Dynamics.Base
 
@@ -48,12 +49,12 @@
 --
 
 -- | Return the maximum.
-maxD :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a
-maxD = liftM2 max
+maxDynamics :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a
+maxDynamics = liftM2 max
 
 -- | Return the minimum.
-minD :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a
-minD = liftM2 min
+minDynamics :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a
+minDynamics = liftM2 min
 
 --
 -- Integrals
@@ -65,10 +66,10 @@
                      integInternal :: IORef (Dynamics Double) }
 
 -- | Create a new integral with the specified initial value.
-newInteg :: Dynamics Double -> Dynamics Integ
+newInteg :: Dynamics Double -> Simulation Integ
 newInteg i = 
-  do r1 <- liftIO $ newIORef $ initD i 
-     r2 <- liftIO $ newIORef $ initD i 
+  do r1 <- liftIO $ newIORef $ initDynamics i 
+     r2 <- liftIO $ newIORef $ initDynamics i 
      let integ = Integ { integInit     = i, 
                          integExternal = r1,
                          integInternal = r2 }
@@ -87,7 +88,7 @@
      m p
 
 -- | Set the derivative for the integral.
-integDiff :: Integ -> Dynamics Double -> Dynamics ()
+integDiff :: Integ -> Dynamics Double -> Simulation ()
 integDiff integ diff =
   do let z = Dynamics $ \p ->
            do y <- readIORef (integExternal integ)
@@ -292,7 +293,7 @@
 -- | Return an integral with the specified derivative and initial value.
 -- If you want to create a loopback then you should use the 'Integ' type 
 -- directly. The 'integ' function is just a wrapper that uses this type.
-integ :: Dynamics Double -> Dynamics Double -> Dynamics (Dynamics Double)
+integ :: Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
 integ diff i =
   do x <- newInteg i
      integDiff x diff
@@ -308,10 +309,10 @@
                    sumInternal :: IORef (Dynamics a) }
 
 -- | Create a new sum with the specified initial value.
-newSum :: (MArray IOUArray a IO, Num a) => Dynamics a -> Dynamics (Sum a)
+newSum :: (MArray IOUArray a IO, Num a) => Dynamics a -> Simulation (Sum a)
 newSum i =   
-  do r1 <- liftIO $ newIORef $ initD i 
-     r2 <- liftIO $ newIORef $ initD i 
+  do r1 <- liftIO $ newIORef $ initDynamics i 
+     r2 <- liftIO $ newIORef $ initDynamics i 
      let sum = Sum { sumInit     = i, 
                      sumExternal = r1,
                      sumInternal = r2 }
@@ -330,7 +331,7 @@
      m p
 
 -- | Set the difference equation for the sum.
-sumDiff :: (MArray IOUArray a IO, Num a) => Sum a -> Dynamics a -> Dynamics ()
+sumDiff :: (MArray IOUArray a IO, Num a) => Sum a -> Dynamics a -> Simulation ()
 sumDiff sum (Dynamics diff) =
   do let z = Dynamics $ \p ->
            case pointIteration p of
diff --git a/Simulation/Aivika/Dynamics/UVar.hs b/Simulation/Aivika/Dynamics/UVar.hs
--- a/Simulation/Aivika/Dynamics/UVar.hs
+++ b/Simulation/Aivika/Dynamics/UVar.hs
@@ -14,17 +14,19 @@
 --
 module Simulation.Aivika.Dynamics.UVar
        (UVar,
-        newUVar,
         uvarQueue,
+        newUVar,
         readUVar,
         writeUVar,
         modifyUVar,
         freezeUVar) where
 
+import Control.Monad
 import Data.Array
 import Data.Array.IO
 import Data.IORef
 
+import Simulation.Aivika.Dynamics.Internal.Simulation
 import Simulation.Aivika.Dynamics.Internal.Dynamics
 import Simulation.Aivika.Dynamics.EventQueue
 
@@ -39,12 +41,12 @@
          uvarYS    :: UV.UVector a}
      
 -- | Create a new variable bound to the specified event queue.
-newUVar :: (MArray IOUArray a IO) => EventQueue -> a -> Dynamics (UVar a)
+newUVar :: (MArray IOUArray a IO) => EventQueue -> a -> Simulation (UVar a)
 newUVar q a =
-  Dynamics $ \p ->
+  Simulation $ \r ->
   do xs <- UV.newVector
      ys <- UV.newVector
-     UV.appendVector xs $ spcStartTime $ pointSpecs p
+     UV.appendVector xs $ spcStartTime $ runSpecs r
      UV.appendVector ys a
      return UVar { uvarQueue = q,
                    uvarRun   = queueRun q,
@@ -121,6 +123,8 @@
               UVar a -> Dynamics (Array Int Double, Array Int a)
 freezeUVar v =
   Dynamics $ \p ->
-  do xs <- UV.freezeVector (uvarXS v)
+  do let Dynamics m = uvarRun v
+     m p
+     xs <- UV.freezeVector (uvarXS v)
      ys <- UV.freezeVector (uvarYS v)
      return (xs, ys)
diff --git a/Simulation/Aivika/Dynamics/Var.hs b/Simulation/Aivika/Dynamics/Var.hs
--- a/Simulation/Aivika/Dynamics/Var.hs
+++ b/Simulation/Aivika/Dynamics/Var.hs
@@ -12,8 +12,8 @@
 --
 module Simulation.Aivika.Dynamics.Var
        (Var,
-        newVar,
         varQueue,
+        newVar,
         readVar,
         writeVar,
         modifyVar,
@@ -23,6 +23,7 @@
 import Data.Array.IO
 import Data.IORef
 
+import Simulation.Aivika.Dynamics.Internal.Simulation
 import Simulation.Aivika.Dynamics.Internal.Dynamics
 import Simulation.Aivika.Dynamics.EventQueue
 
@@ -40,12 +41,12 @@
         varYS    :: V.Vector a}
      
 -- | Create a new variable bound to the specified event queue.
-newVar :: EventQueue -> a -> Dynamics (Var a)
+newVar :: EventQueue -> a -> Simulation (Var a)
 newVar q a =
-  Dynamics $ \p ->
+  Simulation $ \r ->
   do xs <- UV.newVector
      ys <- V.newVector
-     UV.appendVector xs $ spcStartTime $ pointSpecs p
+     UV.appendVector xs $ spcStartTime $ runSpecs r
      V.appendVector ys a
      return Var { varQueue = q,
                   varRun   = queueRun q,
@@ -121,6 +122,8 @@
 freezeVar :: Var a -> Dynamics (Array Int Double, Array Int a)
 freezeVar v =
   Dynamics $ \p ->
-  do xs <- UV.freezeVector (varXS v)
+  do let Dynamics m = varRun v
+     m p
+     xs <- UV.freezeVector (varXS v)
      ys <- V.freezeVector (varYS v)
      return (xs, ys)
diff --git a/Simulation/Aivika/Statistics.hs b/Simulation/Aivika/Statistics.hs
--- a/Simulation/Aivika/Statistics.hs
+++ b/Simulation/Aivika/Statistics.hs
@@ -3,7 +3,7 @@
 
 -- |
 -- Module     : Simulation.Aivika.Statistics
--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- Copyright  : Copyright (c) 2009-2012, David Sorokin <david.sorokin@gmail.com>
 -- License    : BSD3
 -- Maintainer : David Sorokin <david.sorokin@gmail.com>
 -- Stability  : experimental
@@ -16,14 +16,15 @@
         newStatistics,
         addStatistics,
         statisticsData,
-        AnalysisResults(..),
         analyzeData,
+        AnalysisResults(..),
         showResults) where 
 
 import Data.Foldable
 import Data.Array
 import Data.Array.IO
 import Control.Monad
+import Control.Monad.Trans
 import Control.Concurrent.MVar
 
 import Simulation.Aivika.UVector
diff --git a/aivika.cabal b/aivika.cabal
--- a/aivika.cabal
+++ b/aivika.cabal
@@ -1,5 +1,5 @@
 name:            aivika
-version:         0.2
+version:         0.3
 synopsis:        A multi-paradigm simulation library
 description:
     Aivika is a small simulation library that covers many paradigms. 
@@ -7,28 +7,39 @@
     Also it can be applied to the Discrete Event Simulation. It supports 
     the event-oriented, process-oriented and activity-oriented paradigms. 
     Aivika also supports the Agent-based Modeling. Finally, it can be applied 
-    to System Dynamics.
+    to System Dynamics. 
     .
+    It is possible due to using a very general approach when the basic 
+    modeling entity is just a function of simulation time. The paradigms
+    are mainly distinguished by sets of the functions that are used to 
+    model the activities. These sets are small and do not pretend
+    to be comprehensive. Aivika is mostly a proof-of-concept project
+    rather than a big library that knows everything.
+    .
     The library widely uses monads. The dynamic system is represented as 
     a computation in the Dynamics monad. There is also the Process
     monad to represent the discontinuous processes which can suspend
-    at any time and then resume later. Everything else is expressed through 
-    these two monads, including the event handlers, agent handlers and even 
-    integrals.
+    at any time and then resume later. There is also the Simulation monad
+    that represents a simulation run, in which scope the previous 
+    two monads exist. Almost everything is expressed through these monads, 
+    including the event handlers, agent handlers and even integrals 
+    except for the parameters and statistics that already use the IO monad.
     .
     The PDF documentation is available at 
-    <https://github.com/dsorokin/aivika/blob/master/doc/aivika.pdf>
+    <https://github.com/dsorokin/aivika/blob/master/doc/aivika.pdf>.
+    Please note that the documentation is obsolete and it corresponds to 
+    the previous version but it can still be helpful.
     .
 category:        Simulation
 license:         BSD3
 license-file:    LICENSE
-copyright:       (c) 2009-2011. David Sorokin <david.sorokin@gmail.com>
+copyright:       (c) 2009-2012. David Sorokin <david.sorokin@gmail.com>
 author:          David Sorokin
 maintainer:      David Sorokin <david.sorokin@gmail.com>
 homepage:        http://github.com/dsorokin/aivika
 cabal-version:   >= 1.2.0
 build-type:      Simple
-tested-with:     GHC == 7.0.3
+tested-with:     GHC == 7.0.4
 
 extra-source-files:  examples/BassDiffusion.hs
                      examples/ChemicalReaction.hs
@@ -49,20 +60,21 @@
                      Simulation.Aivika.Dynamics.Base
                      Simulation.Aivika.Dynamics.Cont
                      Simulation.Aivika.Dynamics.EventQueue
-                     Simulation.Aivika.Dynamics.Lift
+                     Simulation.Aivika.Dynamics.Parameter
                      Simulation.Aivika.Dynamics.Process
                      Simulation.Aivika.Dynamics.Random
                      Simulation.Aivika.Dynamics.Ref
                      Simulation.Aivika.Dynamics.Resource
+                     Simulation.Aivika.Dynamics.Simulation
                      Simulation.Aivika.Dynamics.SystemDynamics
                      Simulation.Aivika.Dynamics.UVar
                      Simulation.Aivika.Dynamics.Var
-                     Simulation.Aivika.Dynamics.Parameter
+                     Simulation.Aivika.Statistics
                      Simulation.Aivika.PriorityQueue
                      Simulation.Aivika.Queue
-                     Simulation.Aivika.Statistics
 
     other-modules:   Simulation.Aivika.Dynamics.Internal.Dynamics
+                     Simulation.Aivika.Dynamics.Internal.Simulation
                      Simulation.Aivika.Dynamics.Internal.Cont
                      Simulation.Aivika.Dynamics.Internal.Process
                      Simulation.Aivika.Dynamics.Internal.Time
@@ -73,11 +85,12 @@
                      Simulation.Aivika.UVector
                      
     build-depends:   base >= 3 && < 6,
-                     haskell98,
                      mtl >= 1.1.0.2,
                      array >= 0.3.0.0,
-                     containers >= 0.4.0.0
+                     containers >= 0.4.0.0,
+                     random >= 1.0.0.3
 
-    extensions:      FlexibleContexts
+    extensions:      FlexibleContexts,
+                     BangPatterns
                      
     ghc-options:     -O2
diff --git a/examples/BassDiffusion.hs b/examples/BassDiffusion.hs
--- a/examples/BassDiffusion.hs
+++ b/examples/BassDiffusion.hs
@@ -1,10 +1,11 @@
 
-import Random
+import System.Random
 import Data.Array
 import Control.Monad
 import Control.Monad.Trans
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Simulation
 import Simulation.Aivika.Dynamics.EventQueue
 import Simulation.Aivika.Dynamics.Agent
 import Simulation.Aivika.Dynamics.Ref
@@ -34,7 +35,7 @@
                        personPotentialAdopter :: AgentState,
                        personAdopter :: AgentState }
               
-createPerson :: EventQueue -> Dynamics Person              
+createPerson :: EventQueue -> Simulation Person              
 createPerson q =    
   do agent <- newAgent q
      potentialAdopter <- newState agent
@@ -43,14 +44,14 @@
                      personPotentialAdopter = potentialAdopter,
                      personAdopter = adopter }
        
-createPersons :: EventQueue -> Dynamics (Array Int Person)
+createPersons :: EventQueue -> Simulation (Array Int Person)
 createPersons q =
   do list <- forM [1 .. n] $ \i ->
        do p <- createPerson q
           return (i, p)
      return $ array (1, n) list
      
-definePerson :: Person -> Array Int Person -> Ref Int -> Ref Int -> Dynamics ()
+definePerson :: Person -> Array Int Person -> Ref Int -> Ref Int -> Simulation ()
 definePerson p ps potentialAdopters adopters =
   do stateActivation (personPotentialAdopter p) $
        do modifyRef potentialAdopters $ \a -> a + 1
@@ -75,7 +76,7 @@
      stateDeactivation (personAdopter p) $
        modifyRef adopters $ \a -> a - 1
         
-definePersons :: Array Int Person -> Ref Int -> Ref Int -> Dynamics ()
+definePersons :: Array Int Person -> Ref Int -> Ref Int -> Simulation ()
 definePersons ps potentialAdopters adopters =
   forM_ (elems ps) $ \p -> 
   definePerson p ps potentialAdopters adopters
@@ -87,17 +88,20 @@
 activatePersons ps =
   forM_ (elems ps) $ \p -> activatePerson p
 
-model :: Dynamics (Dynamics [Int])
+model :: Simulation [IO [Int]]
 model =
   do q <- newQueue
      potentialAdopters <- newRef q 0
      adopters <- newRef q 0
      ps <- createPersons q
      definePersons ps potentialAdopters adopters
-     activatePersons ps
-     return $ do i1 <- readRef potentialAdopters
-                 i2 <- readRef adopters
-                 return [i1, i2]
+     runDynamicsInStart $
+       activatePersons ps
+     runDynamics $
+       do i1 <- readRef potentialAdopters
+          i2 <- readRef adopters
+          return [i1, i2]
 
-main =
-  printDynamics model specs
+main = 
+  do xs <- runSimulation model specs
+     forM_ xs $ \x -> x >>= print
diff --git a/examples/ChemicalReaction.hs b/examples/ChemicalReaction.hs
--- a/examples/ChemicalReaction.hs
+++ b/examples/ChemicalReaction.hs
@@ -1,5 +1,6 @@
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Simulation
 import Simulation.Aivika.Dynamics.SystemDynamics
 
 specs = Specs { spcStartTime = 0, 
@@ -7,7 +8,7 @@
                 spcDT = 0.01,
                 spcMethod = RungeKutta4 }
 
-model :: Dynamics (Dynamics [Double])
+model :: Simulation [Double]
 model =
   do integA <- newInteg 100
      integB <- newInteg 0
@@ -20,8 +21,6 @@
      integDiff integA (- ka * a)
      integDiff integB (ka * a - kb * b)
      integDiff integC (kb * b)
-     return $ sequence [a, b, c]
+     runDynamicsInFinal $ sequence [a, b, c]
 
-main = 
-  do a <- runDynamics1 model specs
-     print a
+main = runSimulation model specs >>= print
diff --git a/examples/FishBank.hs b/examples/FishBank.hs
--- a/examples/FishBank.hs
+++ b/examples/FishBank.hs
@@ -2,6 +2,7 @@
 import Data.Array
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Simulation
 import Simulation.Aivika.Dynamics.SystemDynamics
 
 specs = Specs { spcStartTime = 0, 
@@ -10,7 +11,7 @@
                 -- spcDT = 0.000005,
                 spcMethod = RungeKutta4 }
 
-model :: Dynamics (Dynamics Double)
+model :: Simulation Double
 model =
   do fishInteg <- newInteg 1000
      shipsInteg <- newInteg 10
@@ -36,23 +37,22 @@
                               (0.6, 5.118), (0.7, 5.247), (0.8, 5.849), 
                               (0.9, 6.151), (10.0, 6.194)]
          density = fish / area
-         fishDeathRate = maxD 0 (fish * deathFraction)
-         fishHatchRate = maxD 0 (fish * hatchFraction)
+         fishDeathRate = maxDynamics 0 (fish * deathFraction)
+         fishHatchRate = maxDynamics 0 (fish * hatchFraction)
          fishPrice = 20
          fractionInvested = 0.2
          hatchFraction = 6
          operatingCost = ships * 250
          profit = revenue - operatingCost
          revenue = totalCatchPerYear * fishPrice
-         shipBuildingRate = maxD 0 (profit * fractionInvested / shipCost)
+         shipBuildingRate = maxDynamics 0 (profit * fractionInvested / shipCost)
          shipCost = 300
-         totalCatchPerYear = maxD 0 (ships * catchPerShip)
+         totalCatchPerYear = maxDynamics 0 (ships * catchPerShip)
      -- derivatives --
      integDiff fishInteg (fishHatchRate - fishDeathRate - totalCatchPerYear)
      integDiff shipsInteg shipBuildingRate
      integDiff totalProfitInteg annualProfit
      -- results --
-     return annualProfit
+     runDynamicsInFinal annualProfit
 
-main = do a <- runDynamics1 model specs
-          print a
+main = runSimulation model specs >>= print
diff --git a/examples/Furnace.hs b/examples/Furnace.hs
--- a/examples/Furnace.hs
+++ b/examples/Furnace.hs
@@ -1,12 +1,12 @@
 
 import Data.Maybe
-import Random
+import System.Random
 import Control.Monad
 import Control.Monad.Trans
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Simulation
 import Simulation.Aivika.Dynamics.Base
-import Simulation.Aivika.Dynamics.Lift
 import Simulation.Aivika.Dynamics.EventQueue
 import Simulation.Aivika.Dynamics.Ref
 import Simulation.Aivika.Dynamics.UVar
@@ -98,7 +98,7 @@
           }
 
 -- | Create a furnace.
-newFurnace :: EventQueue -> Dynamics Furnace
+newFurnace :: EventQueue -> Simulation Furnace
 newFurnace queue =
   do normalGen <- liftIO normalGen
      pits <- sequence [newPit queue | i <- [1..10]]
@@ -133,7 +133,7 @@
                       furnaceUnloadTemps = unloadTemps }
 
 -- | Create a new pit.
-newPit :: EventQueue -> Dynamics Pit
+newPit :: EventQueue -> Simulation Pit
 newPit queue =
   do ingot <- newRef queue Nothing
      h' <- newRef queue 0.0
@@ -255,10 +255,10 @@
      modifyRef (furnaceLoadCount furnace) (+ 1)
   
 -- | Iterate the furnace processing.
-iterateFurnace :: Furnace -> Dynamics (Dynamics ())
+iterateFurnace :: Furnace -> Simulation (Dynamics ())
 iterateFurnace furnace = 
   let pits = furnacePits furnace
-  in iterateD $
+  in iterateDynamics $
      do ready <- ingotsReady furnace
         when ready $ 
           do mapM_ (tryUnloadPit furnace) pits
@@ -303,7 +303,7 @@
   do delay <- liftIO $ exprnd (1.0 / 2.5)
      holdProcess delay
      -- we have got a new ingot
-     liftD $ acceptIngot furnace
+     liftDynamics $ acceptIngot furnace
      -- repeat it again
      processFurnace furnace
 
@@ -338,80 +338,80 @@
     rho x = (x - ex) ^ 2
 
 -- | The simulation model.
-model :: Dynamics (Dynamics ())
+model :: Simulation ()
 model =
   do queue <- newQueue
      furnace <- newFurnace queue
      pid <- newProcessID queue
      
-     initializeFurnace furnace
+     runDynamicsInStart $
+       initializeFurnace furnace
      
      -- get the furnace iterator
      iterator <- iterateFurnace furnace
      
      -- accept input ingots
-     t0 <- starttime
-     runProcess (processFurnace furnace) pid t0
+     runDynamicsInStart $
+       do t0 <- starttime
+          runProcess (processFurnace furnace) pid t0
      
-     let system :: Dynamics ()
-         system = 
-           do iterator   --  iterate in each time point
+     -- run the model in the final time point
+     runDynamicsInFinal $
+       do iterator   --  iterate in each time point
          
-              -- the ingots
-              c0 <- readRef (furnaceTotalCount furnace)
-              c1 <- readRef (furnaceLoadCount furnace)
-              c2 <- readRef (furnaceUnloadCount furnace)
-              c3 <- readRef (furnaceWaitCount furnace)
+          -- the ingots
+          c0 <- readRef (furnaceTotalCount furnace)
+          c1 <- readRef (furnaceLoadCount furnace)
+          c2 <- readRef (furnaceUnloadCount furnace)
+          c3 <- readRef (furnaceWaitCount furnace)
               
-              liftIO $ do
-                putStrLn "The count of ingots:"
-                putStrLn $ "  total  = " ++ show c0
-                putStrLn $ "  loaded = " ++ show c1
-                putStrLn $ "  ready  = " ++ show c2
-                putStrLn $ "  awaited in the queue = " ++ show c3
-                putStrLn ""
+          liftIO $ do
+            putStrLn "The count of ingots:"
+            putStrLn $ "  total  = " ++ show c0
+            putStrLn $ "  loaded = " ++ show c1
+            putStrLn $ "  ready  = " ++ show c2
+            putStrLn $ "  awaited in the queue = " ++ show c3
+            putStrLn ""
          
-              -- the temperature of the ready ingots
-              (n1, e1, d1) <- 
-                fmap stats $ readRef (furnaceUnloadTemps furnace)
+          -- the temperature of the ready ingots
+          (n1, e1, d1) <- 
+            fmap stats $ readRef (furnaceUnloadTemps furnace)
                 
-              liftIO $ do 
-                putStrLn "The temperature of the ready ingots:"
-                putStrLn $ "  average   = " ++ show e1
-                putStrLn $ "  deviation = " ++ show d1
-                putStrLn ""
+          liftIO $ do 
+            putStrLn "The temperature of the ready ingots:"
+            putStrLn $ "  average   = " ++ show e1
+            putStrLn $ "  deviation = " ++ show d1
+            putStrLn ""
                 
-              -- the ingots in pits
-              r2 <- fmap analyzeData $ liftIO $ statisticsData (furnacePitCountStats furnace)
+          -- the ingots in pits
+          r2 <- fmap analyzeData $ liftIO $ statisticsData (furnacePitCountStats furnace)
               
-              liftIO $ do
-                putStrLn "The ingots in pits: "
-                putStrLn $ showResults r2 2 []
-                putStrLn ""
+          liftIO $ do
+            putStrLn "The ingots in pits: "
+            putStrLn $ showResults r2 2 []
+            putStrLn ""
               
-              -- the queue size
-              r3 <- fmap analyzeData $ liftIO $ statisticsData (furnaceQueueCountStats furnace)
+          -- the queue size
+          r3 <- fmap analyzeData $ liftIO $ statisticsData (furnaceQueueCountStats furnace)
      
-              liftIO $ do
-                putStrLn "The queue size: "
-                putStrLn $ showResults r3 2 []
-                putStrLn ""
+          liftIO $ do
+            putStrLn "The queue size: "
+            putStrLn $ showResults r3 2 []
+            putStrLn ""
               
-              -- the mean wait time in the queue
-              t4 <- readRef (furnaceWaitTime furnace) /
-                   fmap (fromInteger . toInteger)
-                   (readRef (furnaceWaitCount furnace))
+          -- the mean wait time in the queue
+          t4 <- readRef (furnaceWaitTime furnace) /
+                fmap (fromInteger . toInteger)
+                (readRef (furnaceWaitCount furnace))
               
-              -- the mean heating time
-              t5 <- readRef (furnaceHeatingTime furnace) /
-                   fmap (fromInteger . toInteger)
-                   (readRef (furnaceUnloadCount furnace))
+          -- the mean heating time
+          t5 <- readRef (furnaceHeatingTime furnace) /
+                fmap (fromInteger . toInteger)
+                (readRef (furnaceUnloadCount furnace))
                     
-              liftIO $ do
-                putStrLn $ "The mean wait time: " ++ show t4
-                putStrLn $ "The mean heating time: " ++ show t5
-         
-     return system
+          liftIO $ do
+            putStrLn $ "The mean wait time: " ++ show t4
+            putStrLn $ "The mean heating time: " ++ show t5
 
 -- | The main program.
-main = runDynamics1 model specs
+main = runSimulation model specs
diff --git a/examples/MachRep1.hs b/examples/MachRep1.hs
--- a/examples/MachRep1.hs
+++ b/examples/MachRep1.hs
@@ -15,12 +15,12 @@
 -- Output is long-run proportion of up time. Should get value of about
 -- 0.66.
 
-import Random
+import System.Random
 import Control.Monad.Trans
 
 import Simulation.Aivika.Dynamics
 import Simulation.Aivika.Dynamics.Base
-import Simulation.Aivika.Dynamics.Lift
+import Simulation.Aivika.Dynamics.Simulation
 import Simulation.Aivika.Dynamics.EventQueue
 import Simulation.Aivika.Dynamics.Ref
 import Simulation.Aivika.Dynamics.Process
@@ -38,7 +38,7 @@
   do x <- getStdRandom random
      return (- log x / lambda)
      
-model :: Dynamics (Dynamics Double)
+model :: Simulation Double
 model =
   do queue <- newQueue
      totalUpTime <- newRef queue 0.0
@@ -48,29 +48,25 @@
      
      let machine :: Process ()
          machine =
-           do startUpTime <- liftD time
+           do startUpTime <- liftDynamics time
               upTime <- liftIO $ exprnd upRate
               holdProcess upTime
-              finishUpTime <- liftD time
-              liftD $ modifyRef totalUpTime
+              finishUpTime <- liftDynamics time
+              liftDynamics $ 
+                modifyRef totalUpTime
                 (+ (finishUpTime - startUpTime))
               repairTime <- liftIO $ exprnd repairRate
               holdProcess repairTime
               machine
          
-     t0 <- starttime
-     
-     runProcess machine pid1 t0
-     runProcess machine pid2 t0
-     
-     let system :: Dynamics Double
-         system =
-           do x <- readRef totalUpTime
-              y <- stoptime
-              return $ x / (2 * y)
+     runDynamicsInStart $
+       do t0 <- starttime
+          runProcess machine pid1 t0
+          runProcess machine pid2 t0
      
-     return system
+     runDynamicsInFinal $
+       do x <- readRef totalUpTime
+          y <- stoptime
+          return $ x / (2 * y)
   
-main =         
-  do a <- runDynamics1 model specs
-     print a
+main = runSimulation model specs >>= print
diff --git a/examples/MachRep1EventDriven.hs b/examples/MachRep1EventDriven.hs
--- a/examples/MachRep1EventDriven.hs
+++ b/examples/MachRep1EventDriven.hs
@@ -15,10 +15,11 @@
 -- Output is long-run proportion of up time. Should get value of about
 -- 0.66.
 
-import Random
+import System.Random
 import Control.Monad.Trans
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Simulation
 import Simulation.Aivika.Dynamics.Base
 import Simulation.Aivika.Dynamics.EventQueue
 import Simulation.Aivika.Dynamics.Ref
@@ -36,7 +37,7 @@
   do x <- getStdRandom random
      return (- log x / lambda)
      
-model :: Dynamics (Dynamics Double)
+model :: Simulation Double
 model =
   do queue <- newQueue
      totalUpTime <- newRef queue 0.0
@@ -62,19 +63,16 @@
               let t = startUpTime + upTime
               enqueue queue t $ machineBroken startUpTime
      
-     t0 <- starttime
-     
-     enqueue queue t0 machineRepaired    -- start the first machine
-     enqueue queue t0 machineRepaired    -- start the second machine
-     
-     let system :: Dynamics Double
-         system =
-           do x <- readRef totalUpTime
-              y <- stoptime
-              return $ x / (2 * y)
-              
-     return system
+     runDynamicsInStart $
+       do t0 <- starttime
+          -- start the first machine
+          enqueue queue t0 machineRepaired
+          -- start the second machine
+          enqueue queue t0 machineRepaired
+          
+     runDynamicsInFinal $
+       do x <- readRef totalUpTime
+          y <- stoptime
+          return $ x / (2 * y)
   
-main =         
-  do a <- runDynamics1 model specs
-     print a
+main = runSimulation model specs >>= print
diff --git a/examples/MachRep1TimeDriven.hs b/examples/MachRep1TimeDriven.hs
--- a/examples/MachRep1TimeDriven.hs
+++ b/examples/MachRep1TimeDriven.hs
@@ -15,10 +15,11 @@
 -- Output is long-run proportion of up time. Should get value of about
 -- 0.66.
 
-import Random
+import System.Random
 import Control.Monad.Trans
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Simulation
 import Simulation.Aivika.Dynamics.Base
 import Simulation.Aivika.Dynamics.EventQueue
 import Simulation.Aivika.Dynamics.Ref
@@ -36,12 +37,12 @@
   do x <- getStdRandom random
      return (- log x / lambda)
      
-model :: Dynamics (Dynamics Double)
+model :: Simulation Double
 model =
   do queue <- newQueue
      totalUpTime <- newRef queue 0.0
      
-     let machine :: Dynamics (Dynamics ())
+     let machine :: Simulation (Dynamics ())
          machine =
            do startUpTime <- newRef queue 0.0 
              
@@ -101,19 +102,14 @@
      m2 <- machine
      
      -- create strictly sequential computations
-     c1 <- iterateD m1
-     c2 <- iterateD m2
-       
-     let system :: Dynamics Double
-         system =
-           do c1    -- involve in the simulation
-              c2    -- involve in the simulation
-              x <- readRef totalUpTime
-              y <- stoptime
-              return $ x / (2 * y)
+     c1 <- iterateDynamics m1
+     c2 <- iterateDynamics m2
      
-     return system
+     runDynamicsInFinal $
+       do c1    -- involve in the simulation
+          c2    -- involve in the simulation
+          x <- readRef totalUpTime
+          y <- stoptime
+          return $ x / (2 * y)
   
-main =         
-  do a <- runDynamics1 model specs
-     print a
+main = runSimulation model specs >>= print
diff --git a/examples/MachRep2.hs b/examples/MachRep2.hs
--- a/examples/MachRep2.hs
+++ b/examples/MachRep2.hs
@@ -17,13 +17,13 @@
 -- that a given machine does not have immediate access to the repairperson 
 -- when the machine breaks down. Output values should be about 0.6 and 0.67. 
 
-import Random
+import System.Random
 import Control.Monad
 import Control.Monad.Trans
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Simulation
 import Simulation.Aivika.Dynamics.Base
-import Simulation.Aivika.Dynamics.Lift
 import Simulation.Aivika.Dynamics.EventQueue
 import Simulation.Aivika.Dynamics.Ref
 import Simulation.Aivika.Dynamics.Resource
@@ -42,7 +42,7 @@
   do x <- getStdRandom random
      return (- log x / lambda)
      
-model :: Dynamics (Dynamics (Double, Double))
+model :: Simulation (Double, Double)
 model =
   do queue <- newQueue
      
@@ -63,18 +63,18 @@
      
      let machine :: Process ()
          machine =
-           do startUpTime <- liftD time
+           do startUpTime <- liftDynamics time
               upTime <- liftIO $ exprnd upRate
               holdProcess upTime
-              finishUpTime <- liftD time
-              liftD $ modifyRef totalUpTime 
+              finishUpTime <- liftDynamics time
+              liftDynamics $ modifyRef totalUpTime 
                 (+ (finishUpTime - startUpTime))
               
               -- check the resource availability
-              liftD $ modifyRef nRep (+ 1)
-              n <- resourceCount repairPerson
+              liftDynamics $ modifyRef nRep (+ 1)
+              n <- liftDynamics $ resourceCount repairPerson
               when (n == 1) $
-                liftD $ modifyRef nImmedRep (+ 1)
+                liftDynamics $ modifyRef nImmedRep (+ 1)
                 
               requestResource repairPerson
               repairTime <- liftIO $ exprnd repairRate
@@ -83,22 +83,17 @@
               
               machine
          
-     t0 <- starttime
-     
-     runProcess machine pid1 t0
-     runProcess machine pid2 t0
-     
-     let system :: Dynamics (Double, Double)
-         system =
-           do x <- readRef totalUpTime
-              y <- stoptime
-              n <- readRef nRep
-              nImmed <- readRef nImmedRep
-              return (x / (2 * y), 
-                      fromIntegral nImmed / fromIntegral n)
-     
-     return system
+     runDynamicsInStart $
+       do t0 <- starttime
+          runProcess machine pid1 t0
+          runProcess machine pid2 t0
+          
+     runDynamicsInFinal $
+       do x <- readRef totalUpTime
+          y <- stoptime
+          n <- readRef nRep
+          nImmed <- readRef nImmedRep
+          return (x / (2 * y), 
+                  fromIntegral nImmed / fromIntegral n)
   
-main =         
-  do a <- runDynamics1 model specs
-     print a
+main = runSimulation model specs >>= print
diff --git a/examples/MachRep3.hs b/examples/MachRep3.hs
--- a/examples/MachRep3.hs
+++ b/examples/MachRep3.hs
@@ -13,13 +13,13 @@
 -- until both machines are down. We find the proportion of up time. It
 -- should come out to about 0.45.
 
-import Random
+import System.Random
 import Control.Monad
 import Control.Monad.Trans
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Simulation
 import Simulation.Aivika.Dynamics.Base
-import Simulation.Aivika.Dynamics.Lift
 import Simulation.Aivika.Dynamics.EventQueue
 import Simulation.Aivika.Dynamics.Ref
 import Simulation.Aivika.Dynamics.Resource
@@ -38,7 +38,7 @@
   do x <- getStdRandom random
      return (- log x / lambda)
      
-model :: Dynamics (Dynamics Double)
+model :: Simulation Double
 model =
   do queue <- newQueue
      
@@ -55,41 +55,38 @@
      
      let machine :: ProcessID -> Process ()
          machine pid =
-           do startUpTime <- liftD time
+           do startUpTime <- liftDynamics time
               upTime <- liftIO $ exprnd upRate
               holdProcess upTime
-              finishUpTime <- liftD time
-              liftD $ modifyRef totalUpTime 
+              finishUpTime <- liftDynamics time
+              liftDynamics $ modifyRef totalUpTime 
                 (+ (finishUpTime - startUpTime))
                 
-              liftD $ modifyRef nUp $ \a -> a - 1
-              nUp' <- liftD $ readRef nUp
+              liftDynamics $ modifyRef nUp $ \a -> a - 1
+              nUp' <- liftDynamics $ readRef nUp
               if nUp' == 1
                 then passivateProcess
-                else do n <- resourceCount repairPerson
-                        when (n == 1) $ reactivateProcess pid
+                else do n <- liftDynamics $ 
+                             resourceCount repairPerson
+                        when (n == 1) $ 
+                          liftDynamics $ reactivateProcess pid
               
               requestResource repairPerson
               repairTime <- liftIO $ exprnd repairRate
               holdProcess repairTime
-              liftD $ modifyRef nUp $ \a -> a + 1
+              liftDynamics $ modifyRef nUp $ \a -> a + 1
               releaseResource repairPerson
               
               machine pid
 
-     t0 <- starttime
-     
-     runProcess (machine pid2) pid1 t0
-     runProcess (machine pid1) pid2 t0
-     
-     let system :: Dynamics Double
-         system =
-           do x <- readRef totalUpTime
-              y <- stoptime
-              return $ x / (2 * y)
+     runDynamicsInStart $
+       do t0 <- starttime
+          runProcess (machine pid2) pid1 t0
+          runProcess (machine pid1) pid2 t0
      
-     return system
+     runDynamicsInFinal $
+       do x <- readRef totalUpTime
+          y <- stoptime
+          return $ x / (2 * y)
   
-main =         
-  do a <- runDynamics1 model specs
-     print a
+main = runSimulation model specs >>= print