diff --git a/src/bin/main.ml b/src/bin/main.ml
index 983c89d..3521d79 100644
--- a/src/bin/main.ml
+++ b/src/bin/main.ml
@@ -7,7 +7,7 @@ open Types
 
 let sample   = ref 10
 let stop     = ref 30.0
-let greedy   = ref false
+let accel    = ref false
 let inplace  = ref false
 let sundials = ref false
 let steps    = ref 1
@@ -26,8 +26,8 @@ let opts = [
   "-sample", Arg.Int (gt0i sample), "n \tSample count (default=10)";
   "-stop", Arg.Float (gt0f stop), "n \tStop time (default=10.0)";
   "-debug", Arg.Set Debug.debug, "\tPrint debug information";
-  "-greedy", Arg.Set greedy, "\tUse greedy simulation";
-  "-sundials", Arg.Set sundials, "\tUse sundials (not compatible with greedy)";
+  "-accelerate", Arg.Set accel, "\tConcatenate continuous functions";
+  "-sundials", Arg.Set sundials, "\tUse sundials (does not support acceleration)";
   "-inplace", Arg.Set inplace, "\tUse imperative solvers";
   "-steps", Arg.Int (gt0i steps), "n \tSplit into [n] steps (default=1)";
 ]
@@ -53,19 +53,18 @@ let st = if !inplace then (module State.InPlaceSimState    : State.SimState)
 
 let sim =
   if !sundials then
-    if !greedy then
-      (Format.eprintf "Sundials does not support greedy simulation\n"; exit 2)
-    else
-      let open StatefulSundials in
-      let c = if !inplace then InPlace.csolve else Functional.csolve in
-      let s = Solver.solver c (d_of_dc z) in
-      let open Sim.LazySim(val st) in run_until_n m s
+    let open StatefulSundials in
+    let c = if !inplace then InPlace.csolve else Functional.csolve in
+    let s = Solver.solver c (d_of_dc z) in
+    let open Sim.Sim(val st) in
+    run_until_n (Output.print !sample (run m s))
   else
     let open StatefulRK45 in
     let c = if !inplace then InPlace.csolve else Functional.csolve in
     let s = Solver.solver_c c z in
-    if !greedy then let open Sim.GreedySim(val st) in run_until_n m s
-               else let open Sim.LazySim(val st) in run_until_n m (d_of_dc s)
+    let open Sim.Sim(val st) in
+    let n = if !accel then accelerate m s else run m (d_of_dc s) in
+    run_until_n (Output.print !sample n)
 
-let () = sim !stop !steps (Output.print !sample)
+let () = sim !stop !steps ignore
 
diff --git a/src/bin/output.ml b/src/bin/output.ml
index a5aceab..c397aa9 100644
--- a/src/bin/output.ml
+++ b/src/bin/output.ml
@@ -1,8 +1,9 @@
 
 open Hsim.Types
+open Hsim.Utils
 open Common
 
-let print_entry t y =
+let print_entry y t =
   let n = Bigarray.Array1.dim y in
   let rec loop i =
     if i = n then ()
@@ -12,17 +13,20 @@ let print_entry t y =
   Format.printf "\n";
   flush stdout
 
-let print samples { h; u } =
+let print_sample samples ({ h; u; _ }, now) =
   let step = h /. (float_of_int samples) in
   let rec loop i =
     if i > samples then ()
-    else if i = samples then print_entry h (u h)
+    else if i = samples then print_entry (u h) (now +. h)
     else let t = float_of_int i *. step in
-    (print_entry t (u t); loop (i+1)) in
-  if h <= 0.0 then begin Debug.print "D: "; print_entry 0.0 (u 0.0) end
+    (print_entry (u t) (now +. t); loop (i+1)) in
+  if h <= 0.0 then begin Debug.print "D: "; print_entry (u 0.0) now end
   else begin Debug.print "C: "; loop 0 end
 
 let print_limits { h; _ } =
   if h <= 0.0 then Format.printf "D: % .10e\n" 0.0
   else Format.printf "C: % .10e\t% .10e\n" 0.0 h
 
+let print samples n =
+  let DNode m = compose n (compose track (map (print_sample samples))) in
+  DNode { m with reset=fun p -> m.reset (p, ((), ())) }
diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml
index 0066293..894aae5 100644
--- a/src/lib/hsim/sim.ml
+++ b/src/lib/hsim/sim.ml
@@ -3,60 +3,70 @@ open Types
 open Solver
 open State
 
-module LazySim (S : SimState) =
+module Sim (S : SimState) =
   struct
     include S
 
-    (** "Lazy" simulation of a model with any solver. *)
+    let step_discrete s step hor fder fzer cget csize zsize jump reset =
+      let ms, ss = get_mstate s, get_sstate s in
+      let i, now, stop = get_input s, get_now s, get_stop s in
+      let o, ms = step ms (i.u now) in
+      let s =
+        let h = hor ms in
+        if h <= 0.0 then set_mstate ms s
+        else if now >= stop then set_idle s
+        else if jump ms then begin
+          let init = cget ms and stop = stop -. now in
+          let fder t = fder ms (Utils.offset i.u now t) in
+          let fzer t = fzer ms (Utils.offset i.u now t) in
+          let ivp = { fder; stop; init; size=csize } in
+          let zc = { init; fzer; size=zsize } in
+          let ss = reset (ivp, zc) ss in
+          let i = { i with h=i.h -. now; u=Utils.offset i.u now } in
+          let mode, stop, now = Continuous, i.h, 0.0 in
+          update ms ss (set_running ~mode ~input:i ~stop ~now s)
+        end else set_running ~mode:Continuous s in
+      Some { h=0.0; c=Discontinuous; u=fun _ -> o }, s
+
+    let step_continuous s step cset fout zset =
+      let ms, ss = get_mstate s, get_sstate s in
+      let i, now, stop = get_input s, get_now s, get_stop s in
+      let (h, f, z), ss = step ss stop in
+      let ms = cset ms (f h) in
+      let fout t = fout ms (i.u (now +. t)) (f (now +. t)) in
+      let s, c = match z with
+      | None ->
+          let s, c = if h >= stop
+            then set_running ~mode:Discrete ~now:h s, Discontinuous
+            else set_running ~now: h s, Continuous in
+          update ms ss s, c
+      | Some z ->
+          let s = set_running ~mode:Discrete ~now:h s in
+          update (zset ms z) ss s, Discontinuous in
+      Some { h=h -. now; u=fout; c }, s
+
+    (** Simulation of a model with any solver. *)
     let run
       (HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
       (DNode solver : ('y, 'yder, 'zin, 'zout) solver)
       : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) lazy_sim
     = let state = get_init model.state solver.state in
+      let step_discrete s =
+        step_discrete s model.step model.horizon model.fder model.fzer
+          model.cget model.csize model.zsize model.jump solver.reset in
 
-      let step s i =
-        let ms, ss = get_mstate s, get_sstate s in
-        match i, is_running s with
-        | Some i, _ ->
+      let step_continuous s =
+        step_continuous s solver.step model.cset model.fout model.zset in
+
+      let step s = function
+        | Some i ->
             let mode, now, stop = Discrete, 0.0, i.h in
-            None, set_running ~mode ~input:i ~now ~stop s
-        | None, false -> None, s
-        | None, true ->
-            let i, now, stop = get_input s, get_now s, get_stop s in
-            match get_mode s with
-            | Discrete ->
-                let o, ms = model.step ms (i.u now) in
-                let s =
-                  let h = model.horizon ms in
-                  if h <= 0.0 then set_mstate ms s
-                  else if now >= stop then set_idle s
-                  else if model.jump ms then begin
-                    let init   = model.cget ms and stop = stop -. now in
-                    let fder t = model.fder ms (Utils.offset i.u now t) in
-                    let fzer t = model.fzer ms (Utils.offset i.u now t) in
-                    let ivp    = { fder; stop; init; size=model.csize } in
-                    let zc     = { init; fzer; size=model.zsize } in
-                    let ss     = solver.reset (ivp, zc) ss in
-                    let i      = { h=i.h -. now; u=Utils.offset i.u now } in
-                    let mode, stop, now = Continuous, i.h, 0.0 in
-                    update ms ss (set_running ~mode ~input:i ~stop ~now s)
-                  end else set_running ~mode:Continuous s in
-                Some { h=0.0; u=fun _ -> o }, s
-            | Continuous ->
-                let (h, f, z), ss = solver.step ss stop in
-                let ms = model.cset ms (f h) in
-                let fout t = model.fout ms (i.u (now +. t)) (f (now +. t)) in
-                let out = { h=h -. now; u=fout } in
-                let s = match z with
-                | None ->
-                    let s = if h >= stop
-                      then set_running ~mode:Discrete ~now:h s
-                      else set_running ~now:h s in
-                    update ms ss s
-                | Some z ->
-                    let s = set_running ~mode:Discrete ~now:h s in
-                    update (model.zset ms z) ss s in
-                Some out, s in
+            step_discrete (set_running ~mode ~input:i ~now ~stop s)
+        | None ->
+            if is_running s then match get_mode s with
+            | Discrete -> step_discrete s
+            | Continuous -> step_continuous s
+            else None, s in
 
       let reset (pm, ps) s =
         let ms = model.reset pm (get_mstate s) in
@@ -65,128 +75,70 @@ module LazySim (S : SimState) =
 
       DNode { state; step; reset }
 
-    (** Run the model on the given input until the end of the input or until the
-        model stops answering. *)
-    let run_on model solver input use =
-      let DNode sim = run model solver in
-      let out = sim.step sim.state (Some input) in
-      let state = match out with None, s -> s | _ -> assert false in
-      let rec loop state =
-        let o, state = sim.step state None in
-        match o with None -> () | Some o -> use o; loop state in
-      loop state
+    let accelerate
+      (HNode  m : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
+      (DNodeC s : ('y, 'yder, 'zin, 'zout) solver_c)
+      : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) lazy_sim
+    = let state = get_init m.state s.state in
+      let step_discrete st =
+        step_discrete st m.step m.horizon m.fder m.fzer m.cget m.csize m.zsize
+          m.jump s.reset in
+      let step_continuous st =
+        step_continuous st s.step m.cset m.fout m.zset in
 
-    (** Run the model on multiple inputs. *)
-    let run_on_n model solver inputs use =
-      let DNode sim = run model solver in
-      ignore @@ List.fold_left (fun state i ->
-        let state = match sim.step state (Some i) with
-        | None, s -> s | _ -> assert false in
-        let rec loop state =
-          let o, state = sim.step state None in
-          match o with None -> state | Some o -> use o; loop state in
-        loop state) sim.state inputs
+      let rec step st = function
+        | Some i ->
+            let mode, now, stop = Discrete, 0.0, i.h in
+            step_discrete (set_running ~mode ~input:i ~now ~stop st)
+        | None ->
+            if is_running st then match get_mode st with
+            | Discrete -> step_discrete st
+            | Continuous ->
+              let o, st = step_continuous st in
+              match o with
+              | None   -> None, st
+              | Some { c=Discontinuous; _ } -> o, st
+              | Some ({ c=Continuous; _ } as o) ->
+                  let o', st = step st None in
+                  match o' with
+                  | None    -> assert false
+                  | Some o' -> Some (Utils.concat [o;o']), st
+            else None, st in
 
-    (** Run the model autonomously until [h], or until the model stops
-        answering. *)
-    let run_until model solver h =
-      run_on model solver { h; u = fun _ -> () }
-
-    (** Run the model autonomously until [length], split in multiple [steps]. *)
-    let run_until_n model solver length steps =
-      let h = length /. float_of_int steps in
-      run_on_n model solver (List.init steps (fun _ -> { h; u=fun _ -> () }))
-  end
-
-module GreedySim (S : SimState) =
-  struct
-    include S
-
-    (** "Greedy" simulation of a model with an appropriate solver. *)
-    let run
-      (HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
-      (DNodeC solver : ('y, 'yder, 'zin, 'zout) solver_c)
-      : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) greedy_sim
-    = let state = get_init model.state solver.state in
-
-      let rec step s i =
-        let ms, ss = get_mstate s, get_sstate s in
-        if not (is_running s) then
-          let mode, now, stop = Discrete, 0.0, i.h in
-          step (set_running ~mode ~input:i ~now ~stop s) i
-        else let now, stop = get_now s, get_stop s in
-        match get_mode s with
-        | Discrete ->
-            let o, ms = model.step ms (i.u now) in
-            let h = model.horizon ms in
-            let rest, s =
-              if h <= 0.0 then step (set_mstate ms s) i
-              else if now >= stop then [], set_idle s
-              else if model.jump ms then
-                let init   = model.cget ms in
-                let fder t = model.fder ms (Utils.offset i.u now t) in
-                let fzer t = model.fzer ms (Utils.offset i.u now t) in
-                let ivp    = { fder; stop = stop -. now; init; size = model.csize } in
-                let zc     = { init; fzer; size = model.zsize } in
-                let ss     = solver.reset (ivp, zc) ss in
-                let i      = { h=i.h -. now; u=Utils.offset i.u now } in
-                let mode, stop, now = Continuous, i.h, 0.0 in
-                step (update ms ss (set_running ~mode ~input:i ~stop ~now s)) i
-              else step (set_running ~mode:Continuous s) i in
-            { h=0.0; u=fun _ -> o }::rest, s
-        | Continuous ->
-            let (h, f, z), ss = solver.step ss stop in
-            let ss = solver.copy ss in
-            let ms = model.cset ms (f h) in
-            let fout t = model.fout ms (i.u (now +. t)) (f (now +. t)) in
-            let out = { h=h -. now; u=fout } in
-            match z with
-            | None ->
-                if h >= stop then
-                  let s = set_running ~mode:Discrete ~now:h s in
-                  let rest, s = step (update ms ss s) i in
-                  out::rest, s
-                else
-                  let s = set_running ~now:h s in
-                  let rest, s = step (update ms ss s) i in
-                  (match rest with
-                   | [] -> [out], s
-                   | f::rest -> Utils.concat [out;f] :: rest, s)
-            | Some z ->
-                let s = set_running ~mode:Discrete ~now:h s in
-                let ms = model.zset ms z in
-                let rest, s = step (update ms ss s) i in
-                out::rest, s in
-
-      let reset (mp, sp) s =
-        let ms = model.reset mp (get_mstate s) in
-        let ss = solver.reset sp (get_sstate s) in
-        update ms ss (set_idle s) in
+      let reset (pm, ps) st =
+        let ms = m.reset pm (get_mstate st) in
+        let ss = s.reset ps (get_sstate st) in
+        update ms ss (set_idle st) in
 
       DNode { state; step; reset }
 
     (** Run the model on the given input until the end of the input or until the
         model stops answering. *)
-    let run_on model solver input use =
-      let DNode sim = run model solver in
-      let o, _ = sim.step sim.state input in
-      List.iter use o
+    let run_on (DNode n) input use =
+      let out = n.step n.state (Some input) in
+      let state = match out with None, s -> s | _ -> assert false in
+      let rec loop state =
+        let o, state = n.step state None in
+        match o with None -> () | Some o -> use o; loop state in
+      loop state
 
     (** Run the model on multiple inputs. *)
-    let run_on_n model solver inputs use =
-      let DNode sim = run model solver in
-      let o, _ = List.fold_left (fun (acc, state) i ->
-        let o, state = sim.step state i in
-        o::acc, state) ([], sim.state) inputs in
-      List.iter use (List.concat (List.rev o))
+    let run_on_n (DNode n) inputs use =
+      ignore @@ List.fold_left (fun state i ->
+        let o, state = n.step state (Some i) in
+        begin match o with None -> () | Some o -> use o end;
+        let rec loop state =
+          let o, state = n.step state None in
+          match o with None -> state | Some o -> use o; loop state in
+        loop state) n.state inputs
 
     (** Run the model autonomously until [h], or until the model stops
         answering. *)
-    let run_until model solver h =
-      run_on model solver { h; u = fun _ -> () }
+    let run_until n h = run_on n { h; c=Discontinuous; u = fun _ -> () }
+
+    (** Run the model autonomously until [h], split in [k] steps. *)
+    let run_until_n n h k =
+      let h = h /. float_of_int k in
+      run_on_n n (List.init k (fun _ -> { h; c=Continuous; u=fun _ -> () }))
 
-    (** Run the model autonomously until [h], split in [n] steps. *)
-    let run_until_n model solver h n =
-      let h = h /. float_of_int n in
-      run_on_n model solver (List.init n (fun _ -> { h; u=fun _ -> () }))
   end
diff --git a/src/lib/hsim/types.ml b/src/lib/hsim/types.ml
index bf9007a..75f7543 100644
--- a/src/lib/hsim/types.ml
+++ b/src/lib/hsim/types.ml
@@ -1,10 +1,12 @@
 
 type time = float
+type continuity = Continuous | Discontinuous
 
 (** Input and output values are functions defined on intervals. *)
 type 'a value =
   { h : time;
-    u : time -> 'a }                   (* Defined on [[0, h]].                *)
+    u : time -> 'a;                    (* Defined on [[0, h]].                *)
+    c : continuity }
 
 (** A time signal is a sequence of possibly absent α-values
     [{ h; u }] where:
diff --git a/src/lib/hsim/utils.ml b/src/lib/hsim/utils.ml
index 0286344..ddeb9bb 100644
--- a/src/lib/hsim/utils.ml
+++ b/src/lib/hsim/utils.ml
@@ -16,7 +16,68 @@ Concatenate functions. [
 let rec concat = function
   | [] -> raise (Invalid_argument "Cannot concatenate an empty value list")
   | [f] -> f
-  | { u; h } :: l ->
-      let { h=hr; u=ur } = concat l in
-      { h=h+.hr; u=fun t -> if t <= h then u t else ur (t -. h) }
+  | { c=Discontinuous; _ } :: _ ->
+      raise (Invalid_argument "Cannot concatenate discontinuous functions")
+  | { u; h; c=Continuous } :: l ->
+      let { h=hr; u=ur; c } = concat l in
+      { c; h=h+.hr; u=fun t -> if t <= h then u t else ur (t -. h) }
 
+let sample { h; u; _ } n =
+  let hs = h /. float_of_int n in
+  let rec step i =
+    if i > n then []
+    else if i = n then [(h, u h)]
+    else let t = float_of_int i *. hs in
+    (t, u t) :: step (i+1) in
+  if h <= 0.0 then [(0.0, u 0.0)] else step 0
+
+(** Compose two nodes together. *)
+let compose (DNode m) (DNode n) =
+  let state = m.state, n.state in
+  let step (ms, ns) i =
+    let v, ms = m.step ms i in
+    let o, ns = n.step ns v in
+    o, (ms, ns) in
+  let reset (ms, ns) (mp, np) =
+    let ms = m.reset ms mp in
+    let ns = n.reset ns np in
+    (ms, ns) in
+  DNode { state; step; reset }
+
+let compose_lazy
+  (DNode m : ('p, 'a, 'b) lazy_sim)
+  (DNode n : ('q, 'b, 'c) lazy_sim)
+= let state = m.state, n.state in
+  let step (ms, ns) = function
+    | Some i ->
+        let v, ms = m.step ms (Some i) in
+        let o, ns = n.step ns v in
+        o, (ms, ns)
+    | None ->
+        let o, ns = n.step ns None in
+        match o with Some o -> Some o, (ms, ns)
+        | None -> let v, ms = m.step ms None in
+        match v with None -> None, (ms, ns)
+        | Some v -> let o, ns = n.step ns (Some v) in
+        o, (ms, ns) in
+  let reset (ms, ns) (mp, np) =
+    let ms = m.reset ms mp in
+    let ns = n.reset ns np in
+    (ms, ns) in
+  DNode { state; step; reset }
+
+(** Track the evolution of a signal in time. *)
+let track : (unit, 'a signal, ('a value * time) option) dnode =
+  let state = 0.0 in
+  let step now = function
+    | None -> None, now
+    | Some i -> Some (i, now), now +. i.h in
+  let reset () _ = 0.0 in
+  DNode { state; step; reset }
+
+(** Apply a function to a signal. *)
+let map f =
+  let state = () in
+  let step () = function None -> None, () | Some v -> Some (f v), () in
+  let reset () () = () in
+  DNode { state; step; reset }