diff --git a/doc/hsim.tex b/doc/hsim.tex
deleted file mode 100644
index b9678f7..0000000
--- a/doc/hsim.tex
+++ /dev/null
@@ -1,145 +0,0 @@
-\documentclass[a4paper]{article}
-\usepackage{fullpage}
-\usepackage{listings}
-\usepackage{xcolor}
-
-\title{A formalization of a simulation engine for hybrid systems}
-\author{}
-\date{}
-
-\setlength{\parindent}{0pt}
-
-\lstdefinelanguage{ml}{
-  basicstyle=\ttfamily,
-  morekeywords=[1]{
-    type, float, option, let, fun, with, if, else, then, in, as, match,
-  },
-  morekeywords=[2]{
-    HNode, DNode, Idle, Running, Discrete, Continuous
-  },
-  keywordstyle=[1]{\color{blue}},
-  keywordstyle=[2]{\color{red}},
-  commentstyle=\itshape,
-  columns=[l]fullflexible,
-  sensitive=true,
-  morecomment=[s]{(*}{*)},
-  keepspaces=true,
-  literate=
-    {'a}{$\alpha$}{1}
-    {'b}{$\beta$}{1}
-    {'p}{$\rho$}{1}
-    {'s}{$\sigma$}{1}
-    {'y}{$y$}{1}
-    {'yder}{$\dot{y}$}{1}
-    {'zin}{$z_{in}$}{1}
-    {'zout}{$z_{out}$}{1}
-    {fun\ }{$\lambda$}{1}
-    {->}{$\to$}{1}
-    {+.}{$+$}{1}
-    {-.}{$-$}{1}
-    {=}{$=$}{1}
-    {>=}{$\geq$}{1}
-    {<=}{$\leq$}{1}
-}
-\lstnewenvironment{ml}{\lstset{language=ml}}{}
-\newcommand{\mlf}[1]{\lstinline[language=ml]{#1}}
-
-
-\begin{document}
-\maketitle
-
-A discrete synchronous function, or node, can be seen as a pair of a step and 
-reset function, which operate on an inner state:
-
-\begin{ml}
-  type ('p, 'a, 'b) dnode =
-    DNode : { s : 's;
-              step  : 's -> 'a -> 'b * 's;
-              reset : 'p -> 's -> 's } -> ('p, 'a, 'b) dnode
-\end{ml}
-
-A hybrid node is quite similar: it has an inner state, a step and a reset 
-function; but the step function is decomposed into multiple distinct elements
-for the purpose of the simulation:
-
-\begin{ml}
-  type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
-    HNode : { s : 's;
-              step : 's -> 'a -> 'b * 's;
-              fder : 's -> 'a -> 'y -> 'yder;
-              fzer : 's -> 'a -> 'y -> 'zout;
-              fout : 's -> 'a -> 'y -> 'b;
-              reset : 'p -> 's -> 's;
-              horizon : 's -> time;
-              jump : 's -> bool;
-              cget : 's -> 'y;
-              cset : 's -> 'y -> 's;
-              zset : 's -> 'zin -> 's;
-            } -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
-\end{ml}
-
-\mlf{step} and \mlf{reset} are on discrete steps, and behave in the same way as 
-for a discrete node. \mlf{fder}, \mlf{fzer} and \mlf{fout} are used during
-integration phases with a solver. \mlf{fder} represents the derivative function
-in an initial value problem. \mlf{fzer} is the zero-crossing function, which
-computes the values of all zero-crossing function we wish to monitor. \mlf{fout}
-computes the output of the system at a particular instant.
-
-\begin{figure}
-   \begin{ml}
-let sim (HNode model) (DNode solver) =
-  let s = { status = Idle; mstate = model.state; sstate = solver.state } in 
-  let step state input =
-    match input, state.status with
-    | Some i, _ ->
-        let status = Running
-          { mode = Continuous; input = i; now = 0.0; stop = i.length } in
-        None, { state with status }
-    | None, Idle -> None, state
-    | None, Running ({ mode = Discrete; _ } as r) ->
-        let o, mstate = model.step state.mstate (r.input.u r.now) in
-        let state =
-          let h = model.horizon mstate in
-          if h <= 0.0 then { state with mstate }
-          else if r.now >= r.stop then s
-          else if model.jump mstate then
-            let y      = model.cget mstate in
-            let fder t = model.fder mstate (offset r.input r.now t) in
-            let fzer t = model.fzer mstate (offset r.input r.now t) in
-            let ivp    = { fder; stop = r.stop -. r.now; init = y } in
-            let zc     = { yc = y; fzer } in
-            let sstate = solver.reset (ivp, zc) state.sstate in
-            let status = Running { r with mode = Continuous } in
-            { status; mstate; sstate }
-          else { state with status = Running { r with mode = Continuous } }
-        in Some { start = r.now; length = 0.0; u = fun _ -> o }, state
-    | None, Running ({ mode = Continuous; _ } as r) ->
-        let (h, f, z), sstate = solver.step state.sstate r.stop in
-        let mstate = model.cset state.mstate (f h) in
-        let now = r.input.start +. h in
-        let state = match z with
-        | None ->
-            let status =
-              if h >= r.stop then Running { r with mode = Discrete; now }
-              else Running { r with now } in
-            { status; mstate; sstate }
-        | Some z ->
-            let status = Running { r with mode = Discrete; now } in
-            { status; mstate = model.zset mstate z; sstate } in
-        let fout t =
-          model.fout mstate (r.input.u (r.now +. t)) (f (r.now +. t)) in
-        let out =
-          { start = r.input.start +. r.now; length = h -. r.now; u = fout } in
-        Some out, state in
-  let reset (m, s) { mstate; sstate; _ } =
-    let mstate = model.reset m mstate in
-    let sstate = solver.reset s sstate in
-    { status = Idle; mstate; sstate } in
-  DNode { state = s; step; reset }
-   \end{ml}
-   \label{fig:ml:sim}
-   \caption{Hybrid System Simulation}
-\end{figure}
-
-
-\end{document}
diff --git a/doc/notes.typ b/doc/notes.typ
index 2e48a1c..6a5cc76 100644
--- a/doc/notes.typ
+++ b/doc/notes.typ
@@ -1,17 +1,36 @@
 #import "@preview/cetz:0.3.2"
 #import "@preview/cetz-plot:0.1.1"
+#import "@preview/finite:0.4.1"
 
 #set page(paper: "a4")
 #set par(justify: true)
+#let breakl(b) = {
+  raw(b.text.replace(regex(" *\\\\\n *"), " "), lang: b.lang, block: b.block)
+}
+
+#show raw: set text(font: "CaskaydiaCove NF")
+#show link: underline
 
 = Notes on hybrid system simulation
 
+== TODO
+
+Multiple different axes of research:
+- abstract over the notion of internal state and of our ability to copy some or
+  part of it
+- make free variables contained in solver closures explicit, and require the
+  ability to copy them for greedy simulation (without the solver itself knowing
+  about them?)
+- think about the notion of input values with overlapping domains
+- think about successive input values carrying information about possible
+  discontinuities
+
 == The problem
 
-The initial problem is _transparent assertions_. We want to simulate assertions 
-with their own solvers, because simulating assertions with the same solver as 
-the parent system changes the final results, when assertions are supposed to be 
-transparent.
+The overarching problem is _transparent assertions_. We want to simulate
+assertions with their own solvers, because simulating assertions with the same
+solver as their parent system changes the final results, when assertions are
+supposed to be transparent.
 
 Assertions are themselves hybrid systems, and can contain their own assertions.
 We thus need a semantics that is both *independent* of the solver, which becomes
@@ -25,73 +44,337 @@ type ('p, 'a, 'b) dnode = DNode :
   { s: 's; step: 's -> 'a -> 'b * 's; reset : 'p -> 's -> 's } ->
   ('p, 'a, 'b) dnode
 ```
-An ODE solver is, in fact, a special case of a discrete node, which is 
+An ODE solver is, in fact, a special case of a discrete node, which is
 initialized with an initial value problem, takes as input a time horizon to
 reach, and returns as output an actual time reached and a dense solution.
 ```ocaml
-type ('y, 'yder) csolver =
-  (('y, 'yder) ivp, time, time * (time -> 'y)) dnode
+type ('y, 'yder) ivp = { init : 'y; fder : time -> 'y -> 'yder; stop : time }
+type ('y, 'yder) csolver = (('y, 'yder) ivp, time, time * (time -> 'y)) dnode
 ```
 Similarly, a zero-crossing solver is initialized with a zero-crossing problem,
 takes as input a horizon and dense function, and returns an actual time reached
 and optional zero-crossing.
 ```ocaml
+type ('y, 'zout) zc = { fzer : time -> 'y -> 'zout; yc : 'y }
 type ('y, 'zin, 'zout) zsolver =
   (('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode
 ```
 A complete solver is then a composition of these two: it is initialized with
 both an `ivp` and `zc`, takes in a horizon, and returns an actual time reached,
 dense solution and optional zero-crossing.
-```ocaml
+
+#breakl(```ocaml
 type ('y, 'yder, 'zin, 'zout) solver =
-  (('y, 'yder) ivp * ('y, 'zout) zc, time, time * (time -> 'y) * 'zin option) dnode
-```
-The construction of a full solver from an ODE and a zero-crossing solver is 
+  (('y, 'yder) ivp * ('y, 'zout) zc, time, time * (time -> 'y) * 'zin option) \
+  dnode
+```)
+
+The construction of a full solver from an ODE and a zero-crossing solver is
 available in the file `src/lib/hsim/solver.ml`.
 
-The simulation of a hybrid system with a solver can itself be seen as a discrete 
+The simulation of a hybrid system with a solver can itself be seen as a discrete
 node, operating on streams of functions defined on successive intervals.
 
-== Output format
+```ocaml
+type ('p, 'a, 'b) sim = ('p, 'a signal, 'b signal) dnode
+```
 
-In "lazy" mode, the output has format ```ocaml 'b signal = 'b value option```.
+== The simulation
 
-In "greedy" mode, the output has format ```ocaml 'b value list```.
+=== Modes and output format
 
-Indeed, we cannot concatenate a discrete step (i.e. a function defined on a 
-singleton interval $[t,t]$) with anything else, as it may be hidden by the
-(inclusive) border of the other function:
+There are two simulation modes, depending on whether the solver provides a way
+to copy its internal state as part of its interface, or more precisely, whether
+the dense solution returned by the solver's step function remains valid after we
+modify the solver's state (through another integration step or a full reset).
+
+In *_lazy_* mode, the output has format ```ocaml 'b signal = 'b value option```.
+The input stream first provides a value for the simulation to run, and must then
+provide as many `None` values as needed for the solver to finish integrating as
+much as possible (that is, until it reaches a zero-crossing or the horizon
+provided by the current input value). This mode is the one used when the solver
+modifies its solution in place (for example, #smallcaps([Sundials])). Each step
+of the simulation runs _one_ step of the solver, and then runs all subsystems as
+much as they need to in order for them to finish using the output of the solver.
+The solver output is then not needed anymore, and we can call the solver again
+to progress with the simulation.
+
+In *_greedy_* mode, the output has format ```ocaml 'b value list```. The input
+stream provides a value for the simulation to run, and the simulation then calls
+the solver as much as it needs to until it reaches the end of the possible
+integration period (zero-crossing or horizon). We concatenate the list of
+results (all functions returned from the solver steps are continuous at their
+edges and defined on more than a single instant, and we can thus create a large,
+piecewise-defined function from all the results up to this point). The solver
+then performs as many discrete steps as needed (i.e. we cascade). These are kept
+in the output list _as is_. We cannot concatenate a discrete step (i.e. a
+function defined on a singleton interval $[t,t]$) with anything else, as it may
+be hidden by the (inclusive) border of the other function. For instance,
 
 #align(center, ```ocaml
 f = concat { start = 0.0; length = 1.0; u = fun t -> t   }
            { start = 1.0; length = 0.0; u = fun t -> 2.0 }
 ```)
 
-#columns(2, [
-  #align(center, cetz.canvas({
-    cetz-plot.plot.plot(
-      size: (1, 1), x-tick-step: 1, y-tick-step: 1, axis-style: "left", {
-        cetz-plot.plot.add(((0,0),(1,1)))
-        cetz-plot.plot.add(((1,2),), mark: "o", mark-size: .03)})}))
-  #colbreak()
-  #linebreak()
+#columns(2, [#{
+  align(center, cetz.canvas({
+    import cetz-plot.plot: *
+    plot(size: (1, 0.5), axis-style: "left", x-tick-step: 1, x-label: "time",
+      y-tick-step: 1, y-label: none, {
+      add(((0,0),(1,1)))
+      add(((1,2),), mark: "o", mark-size: .03)
+    })
+  }))
+  colbreak()
   $ f(t) = cases(
     t   & "if" 0.0 <= t <= 1.0,
     2.0 & "if" t = 1.0 + partial,
     bot & "otherwise",
   ) $
+}])
+
+Here, what is `f.u 1.0`? We cannot differentiate between the multiple
+infinitesimal steps at a given real time `t` using only floating-point numbers
+as input. Therefore, we keep the list of inputs separated, apart from the
+adjacent continuous functions, which we can concatenate together. Until we reach
+the horizon given by the input function, we keep alternating continuous and
+discrete steps. The output is a list of functions, which we can interpret as a
+stream, and the subsystems (the assertions, for instance) can then be called
+repeatedly until they reach the same horizon as the parent system. This has the
+advantage of not needing/having optional values in the input and output streams,
+and the calling system does not need to guess how many steps the simulation will
+need to reach the requested horizon.
+
+=== Steps
+
+A simulation step is either a discrete step or a continuous step.
+
+Discrete steps call the `model.step` function, build a constant output function
+defined on a singleton interval `[now, now]`, and end in one of four possible
+cases:
+- We need to cascade (perform another discrete step)
+- We reached the horizon (reset the solver? idle mode?)
+- We can still integrate, but there has been a state jump (reset solver)
+- We can still integrate, and there has been no state jump (perform a continuous
+  step)
+
+#columns(4, [#align(center, {
+  import cetz: *
+  import cetz-plot.plot: *
+  let p = palette.new(colors: (blue,blue, blue)).with(stroke: true, fill: true)
+  canvas({
+    plot(
+      size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
+      x-ticks: ((1,none),), x-grid: "both", y-tick-step: none, y-label: none,
+      y-min: 1, y-ticks: ((4, none),), y-grid: "both",
+      style: p, {
+      add(((0,2),(.3,2.1),(.7,2.9),(1,3)), line: "spline")
+      add(((1,4),), mark: "o", mark-size: .03)
+      add(((1,5),(1.3,4.7),(1.7,3.2),(2,3)), line: "spline")
+    })
+  })
+  colbreak(); canvas({
+    plot(
+      size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
+      x-ticks: ((3,none),), x-grid: "both", y-tick-step: none, y-label: none,
+      y-min: 1, y-ticks: ((3, none),), y-grid: "both", {
+      add(((0,2),(.5,3.2),(1,5.5),(2,4),(3,3)), line: "spline")
+    })
+  })
+  colbreak(); canvas({
+    plot(
+      size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
+      x-ticks: ((1, none),), x-grid: "both", y-tick-step: none, y-label: none,
+      y-min: 1, y-ticks: ((5, none),), y-grid: "both", {
+      add(((0,2),(0.5,2.3),(1,3)), line: "spline")
+      add(((1,5),(1.5,4.7),(2,4)), line: "spline")
+    })
+  })
+  colbreak(); canvas({
+    plot(
+      size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
+      x-ticks: ((2,none),), x-grid: "both", y-tick-step: none, y-label: none,
+      y-min: 1, y-ticks: ((4, none),), y-grid: "both", {
+      add(((0,2),(0.5,2.3),(1,3),(1.3,5),(2,4),(3,3)), line: "spline")
+    })
+  })
+})])
+
+Continuous steps call the `solver.step` function, parameterized by `model.fder`
+and `model.fzer`, build an output function using `model.fout` defined on the
+interval `[now, now + h]`, and end in one of three possible cases:
+- We found a zero-crossing (perform a discrete step)
+- We found no zero-crossing, but have reached the horizon (perform a discrete
+  step)
+- We found no zero-crossing and have not reached the horizon (perform another
+  continuous step)
+
+#columns(3, [#align(center, {
+  import cetz: *
+  import cetz-plot.plot: *
+  let p = palette.new(colors: (blue,blue, blue)).with(stroke: true, fill: true)
+  canvas({
+    plot(
+      size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
+      x-ticks: ((.7,none),), x-grid: "both", y-tick-step: none, y-label: none,
+      y-min: -2, y-ticks: ((0,$0$),), y-grid: "both", {
+      add(((0,-1),(.5,-.6),(.7,0),(1,2.5),(2,1),(3,.5)), line: "spline")
+    })
+  })
+  colbreak(); canvas({
+    plot(
+      size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
+      x-ticks: ((3,none),), x-grid: "both", y-tick-step: none, y-label: none,
+      y-min: 1, y-ticks: ((3, none),), y-grid: "both", {
+      add(((0,2),(.5,3.2),(1,5.5),(2,4),(3,3)), line: "spline")
+    })
+  })
+  colbreak(); canvas({
+    plot(
+      size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
+      x-ticks: ((2,none),), x-grid: "both", y-tick-step: none, y-label: none,
+      y-min: 1, y-ticks: ((4, none),), y-grid: "both", {
+      add(((0,2),(0.5,2.3),(1,3),(1.3,5),(2,4),(3,3)), line: "spline")
+    })
+  })
+})])
+
+=== Resets
+
+==== Solver reset
+
+When do we need to reset the solver itself?
+
+If we receive a new input that is not continuous w.r.t. the previous input, we
+should reset the solver before starting the integration. This implies that an
+input value is tagged with information relating it to the previous input value.
+If this is not the case, we can reset whenever we receive a new input value,
+although this is less satisfactory.
+
+If a discrete step results in a state jump, we should reset the solver before
+continuing the integration.
+
+==== Simulation reset
+
+Two possible options for the simulation reset:
+- ask for reset parameters for both the model and solver:
+  ```ocaml
+  'p * (('y, 'yder) ivp *  ('y, 'zin, 'zout) zc
+  ```
+- ask for only what is needed for the reset of the model, then build a new
+  initial value problem and zero-crossing problem using the new model state and
+  reset the solver using them. We can build `fder` and `fzer` quite easily. We
+  still need to build `init`, `yc` and `time`.
+
+  `init` and `yc` need to use `model.cget` in order to obtain a continuous state
+  element. However, the usual simulation process first starts by executing a
+  discrete step to initialize the model state. Does `model.cget` make sense
+  before this first step has been executed?
+  #figure(finite.automaton(
+    (D: (D: "cascade", C: "no cascade"),
+     C: (C: "no zero-crossing", D: "zero-crossing")),
+    initial: "D", final: (), layout: finite.layout.linear.with(spacing: 3)
+  )) <automaton>
+  Tagged values that indicate whether they are discontinuous w.r.t. the previous
+  input value do not help us here. Indeed, we could consider that the first
+  value in the stream is always considered discontinuous, which would force a
+  reset as we receive an input value. However, the initialization of the solver
+  still depends on the initialization of the model, which (might) need a
+  discrete step, and we are back to square one.
+
+  Additionally, what stop time do we use for the initial state?
+
+  If possible, this would have the advantage of making the reset function only
+  dependent on the model, since the solver is updated to match the model's
+  simulation requirements.
+  ```ocaml
+  val sim : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode ->
+            ('y, 'yder, 'zin, 'zout) solver ->
+            ('p, 'a, 'b) sim
+  ```
+
+== Mathematical model
+
+#link("https://zelus.di.ens.fr/cc2015/fullpaper.pdf")[[CC'15]] defines the
+simulation loop as follows:
+
+#block(fill: rgb(230, 230, 230), stroke: black, inset: 5pt, [
+  The simulation loop of a hybrid system is a _synchronous function_ defining
+  four streams $t(n)$, $l x(n)$, $y(n)$ and $z(n)$ with $n in bb(N)$.
+  $t(n) in bb(R)$ is the increasing sequence of instants at which the solver
+  stops. $l x(n)$ is the value at time $t(n)$ of the _continuous state
+  variables_, that is, of all variables defined by their derivatives in the
+  original model. $y(n)$ is the value at time $t(n)$ of the _discrete state_.
+  $z(n)$ indicates any _zero-crossings_ at instant $t(n)$ on signals monitored
+  by the solver, that is, any signals that become equal to or pass through zero.
+
+  The synchronous function has two modes: the discrete mode ($D$) contains all
+  computations that may change the discrete state or that have side effects. The
+  continuous mode ($C$) is where ODEs are solved. The two modes alternate
+  according to the execution scheme summarized in @automaton.
+
+  In the _continuous mode_, the solver computes an approximation of the solution
+  of the ODEs and monitors a set of expressions for zero-crossings. The solver
+  is abstracted by a function $s o l v e(f)(g)$ parameterized by $f$ and $g$
+  where:
+  - $x'(tau)=f(y(n),tau,x(tau))$ defines the derivatives of continuous state
+    variables $x$ at instant $tau in bb(R)$;
+  - $u p z(tau)=g(y(n),tau,x(tau))$ defines the current values of a set of
+    zero-crossing signals $u p z$, indexed by $i in {1, ..., k}$.
+  The continuous mode $C$ computes four sequences $s$, $l x$, $z$ and $t$ such
+  that $ (l x, z, t, s)(n+1) = s o l v e(f)(g)(s,y,l x,t,s t e p)(n) $ where
+  - $s(n)$ is the internal state of the solver at instant $t(n) in bb(R)$.
+    Calling $s o l v e(f)(g)$ updates the state to $s(n+1)$.
+  - $x$ is an approximation of a solution of the ODE
+    $ x'(tau)=f(y(n),tau,x(tau)) $
+    It is parameterized by the current discrete state $y(n)$ and initialized at
+    instant $t(n)$ with the value of $l x(n)$, that is, $x(t(n)) = l x(n)$.
+  - $l x(n+1)$ is bounded by the horizon $t(n) + s t e p(n)$ that the solver has
+    been asked to reach, that is
+    $ t(n) <= t(n+1) <= t(n) + s t e p(n) $
+  - $z(n+1)$ signals any zero-crossings detected at time $t(n+1)$. An event
+    occurs with a transition to the discrete mode $D$ when horizon
+    $t(n) + s t e p(n)$ is reached, or when at least one of the zero-crossing
+    signals $u p z(i)$, for $i in {1, ..., k}$ crosses zero, which is indicated
+    by a true value for the corresponding boolean output $z(n+1)(i)$. $
+    e v e n t = z(n+1)(0) or ... or z(n+1)(k) or (t(n+1) = t(n) + s t e p(n)) $
+  If the solver raises an error (for example, a division by zero or an inability
+  to find a solution), we consider that the simulation fails.
+
+  All discrete changes occur in the _discrete mode_ $D$. It is entered when an
+  event is raised during integration. During a discrete phase, the function
+  $n e x t$ defines $y$, $l x$, $s t e p$, $e n c o r e$, $z$ and $t$:
+  $ (y, l x, s t e p, e n c o r e)(n+1) & = n e x t(y, l x, z, t)(n) \
+                                 z(n+1) & = f a l s e \
+                                 t(n+1) & = t(n) $
+  where
+  - $y(n+1)$ is the new discrete state; outside of mode $D$, $y(n+1) = y(n)$.
+  - $l x(n+1)$ is the new continuous state, which may be changed directly in the
+    discrete mode.
+  - $s t e p(n+1)$ is the new step size.
+  - $e n c o r e(n+1)$ is true if an additional discrete step must be performed.
+    Function $n e x t$ can decide to trigger instantaneously another discrete
+    event, causing an _event cascade_.
+  - $t(n)$ (the simulation time) is unchanged during discrete phases.
+  The initial values for $y(0)$, $l x(0)$ and $s(0)$ are given by an
+  initialization function $i n i t$. Finally, $s o l v e(f)(g)$ may decide to
+  reset its internal state if the continuous state changes. If
+  $i n i t \_ s o l v e$ initializes the solver state, we have
+  $ r e i n i t & = (l x(n+1) != l x(n)) \
+         s(n+1) & = "if" r e i n i t "then" i n i t \_ s o l v e(l x(n+1),s(n))
+                    "else" s(n) $
+  Taken together, the definitions from both modes give a synchronous
+  interpretation of the simulation loop as a stream function that computes the
+  sequences $l x$, $y$ and $t$ at instant $n+1$ according to their values at
+  instant $n$ and an internal state. Writing $s o l v e(f)(g)$ abstracts from
+  the actual choice of integration method and zero-crossing detection algorithm.
+  A more detailed description of $s o l v e(f)(g)$ would be possible (for
+  example, an automaton with two states: one that integrates, and one that
+  detects zero-crossings) but with no influence on the code generation problem
+  which must be independent of such simulation details.
+
+  Given a program written in a high-level language, we must produce the
+  functions $i n i t$, $f$, $g$, and $n e x t$. In practice, they are
+  implemented in an imperative language like C. Code generation for hybrid
+  models has much in common with code generation for synchronous languages.
 ])
-
-How do we represent infinitesimal steps in floating-point numbers ?
-
-A possible representation of the above could be
-#align(center, ```ocaml
-f = [{ start = 0.0; length = 1.0; u = fun t -> t   };
-     { start = 1.0; length = 0.0; u = fun _ -> 2.0 }]
-```)
-but returning a list is problematic, however, as it cannot be passed easily as 
-input to the underlying models. Unless the input can itself be considered as a 
-list, in which case the simulation simply operates over the whole list, 
-resetting the solver as needed. This has the advantage of not having to wait 
-several steps, providing nothing, until the solver is done integrating the 
-input, before we can provide it with another input.
diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml
index c7222c7..6c05702 100644
--- a/src/lib/hsim/sim.ml
+++ b/src/lib/hsim/sim.ml
@@ -6,19 +6,28 @@ open State
 let offset (input : 'a value) (now : time) : time -> 'a =
   fun t -> input.u ((now -. input.start) +. t)
 
+let rec compose = function
+  | [] -> assert false
+  | [f] -> f
+  | { start=sl; u=ul; _ } :: l ->
+      let { start=sr; length=lr; u=ur } = compose l in
+      let length = sr +. lr -. sl in
+      { start=sl; length; u=fun t -> if t <= sr then ur t else ul t }
+
 module LazySim (S : SimState) =
   struct
 
-    (* TODO: figure out where we initialize the solvers; the initialization
-             function already supposes an initialized solver state, but could we 
-             parameterize [LazySim] with a solver state module that provides its 
-             own initialization function ? *)
+    (* TODO: figure out where we initialize the solvers; the initialization and
+             reset functions already suppose an initialized solver state, but
+             could we parameterize simulation with a solver state module that
+             provides its own initialization function ? *)
 
     let sim
       (HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
       (DNode solver : ('y, 'yder, 'zin, 'zout) solver)
-      : ('p, 'a, 'b) sim
-    = let s = S.init ~mstate:model.state ~sstate:solver.state in
+      : ('p, 'a, 'b) lazy_sim
+    = let state = S.init ~mstate:model.state ~sstate:solver.state in
+
       let step state input =
         let mstate = S.mstate state and sstate = S.sstate state
         and status = S.status state in
@@ -39,8 +48,9 @@ module LazySim (S : SimState) =
                   let h = model.horizon mstate in
                   if h <= 0.0 then S.update ~mstate state
                   else if now >= stop then
-                    (* TODO: equivalent of [s] (initial state of model and
-                             solvers) ? *)
+                    (* TODO: Build an initial state with the initial states for
+                             both the solver and model - an equivalent of [s] in
+                             the original version. *)
                     raise Common.Utils.TODO
                   else if model.jump mstate then
                     let y      = model.cget mstate in
@@ -75,77 +85,97 @@ module LazySim (S : SimState) =
                     let mstate = model.zset mstate z in
                     S.update ~status ~mstate ~sstate state in
                 Some out, state in
-      let reset m s =
-        let mstate = model.reset m (S.mstate s) in
+
+      let reset p s =
+        (* TODO: does [model.cget mstate] make sense before the first discrete
+                 step - can we use [mstate] to reinitialize [sstate] ? *)
+        let mstate = model.reset p (S.mstate s) in
         let y = model.cget mstate in
+        (* TODO: what initial stop time do we use ? *)
         let stop = raise Common.Utils.TODO in
         let ivp = { fder = model.fder mstate; stop; init = y } in
         let zc = { fzer = model.fzer mstate; yc = y } in
         let sstate = solver.reset (ivp, zc) (S.sstate s) in
         let status = S.idle (S.status s) in
         S.update ~status ~mstate ~sstate s in
-      DNode { state = s; step; reset }
+
+      DNode { state; step; reset }
 
   end
 
 module GreedySim (S : SimState) =
   struct
 
-    (* TODO: greedy simulation *)
+    (* TODO: greedy simulation: call the solvers and the subsystems as often as
+             needed until we reach the horizon. *)
 
     let sim
       (HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
       (DNode solver : ('y, 'yder, 'zin, 'zout) solver)
-      : ('p, 'a, 'b) sim
+      : ('p, 'a, 'b) greedy_sim
     = let state = S.init ~mstate:model.state ~sstate:solver.state in
+
       let rec step state input =
         let status = S.status state and mstate = S.mstate state
         and sstate = S.sstate state in
-        match input, S.is_running state with
-        | Some input, _ ->
-            let mode = Discrete and now = 0.0 and stop = input.length in
-            let status = S.running ~mode ~input ~now ~stop (S.status state) in
-            let state = S.update ~status state in
-            step state None
-        | None, false -> None, state
-        | None, true ->
-            let input = S.input state and now = S.now state
-            and stop = S.stop state in
-            match S.mode state with
-            | Discrete ->
-                let o, mstate = model.step mstate (input.u now) in
-                let state =
-                  let h = model.horizon mstate in
-                  if h <= 0.0 then S.update ~mstate state
-                  else if now >= stop then raise Common.Utils.TODO
-                  else if model.jump mstate then
-                    let y = model.cget mstate in
-                    let fder t = model.fder mstate (offset input now t) in
-                    let fzer t = model.fzer mstate (offset input now t) in
-                    let ivp = { fder; stop = stop -. now; init = y } in
-                    let zc = { yc = y; fzer } in
-                    let sstate = solver.reset (ivp, zc) sstate in
-                    let status = S.running ~mode:Continuous status in
-                    S.update ~status ~mstate ~sstate state
-                  else
-                    let status = S.running ~mode:Continuous status in
-                    S.update ~status state in
-                let start = input.start +. now in
-                Some { start; length = 0.0; u = fun _ -> o }, state
-            | Continuous ->
-                let (h, f, z), sstate = solver.step sstate stop in
-                let mstate = model.cset mstate (f h) in
-                let h' = input.start +. h in
-                let fout t =
-                  model.fout mstate (input.u (now +. t)) (f (now +. t)) in
-                let out =
-                  { start = input.start +. now; length = h -. now; u = fout } in
-                match z with
-                | None ->
-                    let status =
-                      if h >= stop then S.running ~mode:Discrete ~now:h' status
-      in
+        if not (S.is_running state) then
+          let mode = Discrete and now = 0.0 and stop = input.length in
+          let status = S.running ~mode ~input ~now ~stop (S.status state) in
+          let state = S.update ~status state in
+          step state input
+        else let now = S.now state and stop = S.stop state in
+        match S.mode state with
+        | Discrete ->
+            let o, mstate = model.step mstate (input.u now) in
+            let h = model.horizon mstate in
+            let rest, state =
+              if h <= 0.0 then step (S.update ~mstate state) input
+              else if now >= stop then [], state
+              else if model.jump mstate then
+                let y = model.cget mstate in
+                let fder t = model.fder mstate (offset input now t) in
+                let fzer t = model.fzer mstate (offset input now t) in
+                let ivp = { fder; stop = stop -. now; init = y } in
+                let zc = { yc = y; fzer } in
+                let sstate = solver.reset (ivp, zc) sstate in
+                let status = S.running ~mode:Continuous status in
+                step (S.update ~status ~mstate ~sstate state) input
+              else
+                let status = S.running ~mode:Continuous status in
+                step (S.update ~status state) input in
+            let start = input.start +. now in
+            { start; length = 0.0; u = fun _ -> o }::rest, state
+        | Continuous ->
+            let (h, f, z), _sstate = solver.step sstate stop in
+            let mstate = model.cset mstate (f h) in
+            let h' = input.start +. h in
+            let fout t =
+              model.fout mstate (input.u (now +. t)) (f (now +. t)) in
+            let out =
+              { start = input.start +. now; length = h -. now; u = fout } in
+            match z with
+            | None ->
+                if h >= stop then
+                  let status = S.running ~mode:Discrete ~now:h' status in
+                  let rest, state =
+                    step (S.update ~status ~mstate ~sstate state) input in
+                  out::rest, state
+                else
+                  let status = S.running ~now:h' status in
+                  let rest, state =
+                    step (S.update ~status ~mstate ~sstate state) input in
+                  (match rest with
+                   | [] -> [out], state
+                   | f::rest -> compose [out;f] :: rest, state)
+            | Some z ->
+                let status = S.running ~mode:Discrete ~now:h' status in
+                let mstate = model.zset mstate z in
+                let rest, state =
+                  step (S.update ~status ~mstate ~sstate state) input in
+                out::rest, state in
+
       let reset = assert false in
+
       DNode { state; step; reset }
 
   end
diff --git a/src/lib/hsim/state.ml b/src/lib/hsim/state.ml
index 2827b73..28a8856 100644
--- a/src/lib/hsim/state.ml
+++ b/src/lib/hsim/state.ml
@@ -1,20 +1,6 @@
 
 open Types
 
-(* Model *)
-
-module type ModelState =
-  sig
-    type state
-  end
-
-(* Solvers *)
-
-module type SolverState =
-  sig
-    type state
-  end
-
 (* Simulation *)
 
 (** Step mode: discrete or continuous *)
@@ -52,7 +38,7 @@ module type SimState =
     val idle : 'a status -> 'a status
 
     (** Update the status to running. *)
-    val running : 
+    val running :
       ?mode:mode -> ?input:'a value ->
         ?now:time -> ?stop:time -> 'a status -> 'a status
 
@@ -81,7 +67,7 @@ module FunctionalSimState : SimState =
 
     (** Simulation status:
       - [Idle]:    Waiting for input, no activity;
-      - [Running]: Currently integrating: step [mode], current [input], [now] 
+      - [Running]: Currently integrating: step [mode], current [input], [now]
                    timestamp, and [stop] time. *)
     type 'a status =
       | Idle    : 'a status
@@ -92,7 +78,7 @@ module FunctionalSimState : SimState =
             stop  : time;                    (** How long until we stop.      *)
           } -> 'a status
 
-    (** Internal state of the simulation node: model state, solver state and 
+    (** Internal state of the simulation node: model state, solver state and
         current simulation status. *)
     type ('a, 'ms, 'ss) state =
       { status : 'a status;                  (** Current simulation status.   *)
@@ -137,7 +123,7 @@ module FunctionalSimState : SimState =
       match s.status with Running r -> r.input | Idle -> raise Not_running
     let now s =
       match s.status with Running r -> r.now   | Idle -> raise Not_running
-    let stop s = 
+    let stop s =
       match s.status with Running r -> r.stop  | Idle -> raise Not_running
 
     let init ~mstate ~sstate = { status = Idle; mstate; sstate }
@@ -198,7 +184,7 @@ module InPlaceSimState : SimState =
       match s.status with Running r -> r.input | Idle -> raise Not_running
     let now s =
       match s.status with Running r -> r.now   | Idle -> raise Not_running
-    let stop s = 
+    let stop s =
       match s.status with Running r -> r.stop  | Idle -> raise Not_running
 
     let init ~mstate ~sstate = { status = Idle; mstate; sstate }
diff --git a/src/lib/hsim/types.mli b/src/lib/hsim/types.mli
index b181874..64e014e 100644
--- a/src/lib/hsim/types.mli
+++ b/src/lib/hsim/types.mli
@@ -81,5 +81,10 @@ type ('y, 'yder, 'zin, 'zout) solver =
 
 (** The simulation of a hybrid system is a synchronous function on streams of
     functions. *)
-type ('p, 'a, 'b) sim =
+type ('p, 'a, 'b) lazy_sim =
   ('p, 'a signal, 'b signal) dnode
+
+(** Greedy simulation takes in an input and computes as many solver and
+    subsystem steps as needed to reach the input's horizon. *)
+type ('p, 'a, 'b) greedy_sim =
+  ('p, 'a value, 'b value list) dnode