diff --git a/doc/img/ball.png b/doc/img/ball.png
new file mode 100644
index 0000000..2d578f4
Binary files /dev/null and b/doc/img/ball.png differ
diff --git a/doc/notes.typ b/doc/notes.typ
index 3bc8777..ecc9b8f 100644
--- a/doc/notes.typ
+++ b/doc/notes.typ
@@ -129,13 +129,11 @@ cases:
 #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, {
+      y-min: 1, y-ticks: ((4, none),), y-grid: "both", {
       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")
@@ -180,7 +178,6 @@ interval `[now, now + h]`, and end in one of three possible cases:
 #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,
@@ -319,6 +316,71 @@ 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.
 
+=== Composing simulations
+
+The composition of two simulations must maintain the following invariant: the
+output is only ever absent if the solvers have finished integrating up to the
+requested horizon. This is because we need to keep providing absent values
+(`None`) to the first simulation until it reaches the requested horizon,
+otherwise we reset beforehand, and the only information we have on whether the
+simulation has reached its horizon is the nature of its output. Naïve
+composition breaks this:
+
+#grid(columns: (1fr, 1fr), align: center + horizon, cetz.canvas({
+  import cetz.draw: *
+  rect((0, 0), (1, 1))
+  rect((2, 0), (3, 1))
+  content((0.5, 0.5), $M$)
+  content((2.5, 0.5), $N$)
+  line((-0.5, 0.5), (0, 0.5), mark: (end: "straight"))
+  line((1, 0.5), (2, 0.5), mark: (end: "straight"))
+  line((3, 0.5), (3.5, 0.5), mark: (end: "straight"))
+}),
+$ M_i & : & a_1 & | & bot                     & |
+\ M_o & : & b_1 & | & b_2                     & |
+\ N_i & : & b_1 & | & b_2                     & |
+\ N_o & : & c_1 & | & #text(fill: red, $c_2$) & |
+$)
+
+Instead, we need to make sure we only provide $N$ with a new value when
+needed. This is done through the following composition:
+#grid(columns: (1fr, 1fr), align: center + horizon, cetz.canvas({
+  import cetz.draw: *
+  rect((0, 0), (1, 1))
+  rect((2, 0), (3, 1))
+  content((0.5, 0.5), $M$)
+  content((2.5, 0.5), $N$)
+  line((-0.5, 0.5), (0, 0.5), mark: (end: "straight"))
+  line((1, 0.5), (2, 0.5), mark: (end: "straight"))
+  line((3, 0.5), (3.5, 0.5), mark: (end: "straight"))
+}),
+$ I   &: & a_1 & | & bot & | & bot &   &     & | & bot & | & bot &   &     & |
+\ M_i &: & a_1 & | &     & | &     & | & bot & | &     & | &     & | & bot & |
+\ M_o &: & b_1 & | &     & | &     & | & b_2 & | &     & | &     & | & bot & |
+\ N_i &: & b_1 & | & bot & | & bot & | & b_2 & | & bot & | & bot & | & bot & |
+\ N_o &: & c_1 & | & c_2 & | & bot & | & c_3 & | & c_4 & | & bot & | & bot & |
+\ O   &: & c_1 & | & c_2 & | &     &   & c_3 & | & c_4 & | &     &   & bot & |
+$)
+
+```ocaml
+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)
+```
+
 === Resets
 
 ==== Solver reset
@@ -411,6 +473,67 @@ Two possible options for the simulation reset:
   makes this impossible. We thus need reset parameters for both the model and
   solver.
 
+=== Assertions
+
+Assertions are themselves hybrid systems, with their own internal state, `fder`,
+`fzer`, `fout` and `reset` functions, etc.
+A hybrid model with assertions is then a recursive datatype
+```ocaml
+type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a =
+  HNodeA : {
+    body : ('s, 'p,' a, 'b, 'y, 'yder, 'zin, 'zout) hrec;
+    assertions : ('p, 's, float, 'y, 'yder, 'zin, 'zout) hnode_a list;
+  }
+```
+
+Assertions in sub-models need to be "lifted" to the parent model in order to be
+visible for the simulation. Since the inner state of sub-models is contained in
+the state of the parent model, these assertions can be composed with projectors
+on the inner state of the parent model in order to provide them with the correct
+state.
+
+What does it mean for an assertion to run in continuous time ?
+There are two possible approaches:
+- Sampling: take "snapshots" of the state at various timestamps and check the
+  assertions at these snapshots.
+  - Problems: expensive memory usage, and requires state copies for the model.
+    Can also impact stop times for the parent model, which defeats the point of
+    transparent assertions.
+- Real-time: consider the solution provided by the solver as a continuous
+  function and use it as input to monitor for zero-crossings
+  - Problems: limits the expressions we can use in continuous assertions to
+    continuous functions; boolean expressions are inherently discontinuous. We
+    can model inequalities using zero-crossing constructs like `up`
+    (for instance, `a <= b` $-->$ `up(a -. b)`).
+
+Building a continuous function of the entire state is possible. Since the entire
+state can be considered as the sum of a discrete and a continuous part, and that
+the discrete part is considered piecewise-constant, we can build a function
+`time -> 's` using the model's `cset` function (which updates the state's
+continuous part) as follows:
+#grid(columns: (1fr, 1fr), align: center + horizon, ```ocaml
+let f_st (ms : 's) (f : time -> 'y) =
+  fun t -> cset ms (f t)
+```, cetz.canvas({
+  import cetz-plot.plot: *
+  plot(
+    axis-style: "left", x-label: none, x-tick-step: none,
+    y-label: none, y-tick-step: none, y-min: 0.1, {
+    add(((0, 1), (1, 1)), style: (stroke: blue), label: $D$)
+    add(((1, 2), (2.5, 2)), style: (stroke: blue))
+    add(domain: (0, 1), t => - calc.sin(t + 1) + 1.5, style: (stroke: red))
+    add(domain: (1, 2.5), t => calc.sin(t) + 2, style: (stroke: red), label: $C$)
+  })
+}))
+
+If `cset` is pure (that is, it returns a copy of the state, rather than modify
+it in-place), this function poses no problem. If `cset` is not pure and modifies
+the state in-place, this is still useable as long as `f_state` is never used in
+parallel: since `cset` only modifies the continuous part of the state, we can
+always recover the original state by calling `f_state` with the appropriate
+horizon. However, the model's continuous state needs to be re-updated to the
+reached horizon after all assertions/submodels have finished running.
+
 == Mathematical model
 
 #link("https://zelus.di.ens.fr/cc2015/fullpaper.pdf")[[CC'15]] defines the
@@ -497,6 +620,70 @@ simulation loop as follows:
   models has much in common with code generation for synchronous languages.
 ])
 
+In [Branicky, 1995b], dynamical systems are defined as follows:
+
+#block(fill: rgb(230, 230, 230), stroke: black, inset: 5pt, [
+  A continuous (resp. discrete) dynamical system defined on the topological
+  space $X$ over the semigroup $bb(R)^+$ (resp. $bb(Z)^+$) is a function
+  $f : X times bb(R)^+ -> X$ (resp. $f : X times bb(Z)^+ -> X$) with the
+  following three properties:
+  1. initial condition: $f(p, 0) = p$,
+  2. continuity on both arguments,
+  3. semigroup property: $f(f(p, t_1), t_2) = f(p, t_1 + t_2)$
+  for any point $p in X$ and any $t_1, t_2 in bb(R)^+$ (resp $bb(Z)^+$).
+
+  Technically, such functions are referred to as semidynamical systems, with the
+  term dynamical system reserved for those which the semigroups $bb(R)^+$,
+  $bb(Z)^+$ above may be replaced by the groups $bb(R)$, $bb(Z)$. However, the
+  more "popular" notion of dynamical system in math and engineering - and the
+  one used here - requires only the semigroup property. Thus, the term
+  _reversible_ dynamical system is used when it is necessary to distinguish the
+  group from semigroup case.
+
+  The shorthand $[X, S, f]$ denotes the dynamical system $f$ defined on $X$ over
+  the semigroup $S$; $X$ is referred to as its _state space_ and points in $X$
+  are called _states_. If a dynamical system is defined on a subset of $X$, we
+  say it is a dynamical system _in_ $X$.
+
+  For every fixed value of the parameter $s$, the function $f(dot, s)$ defines
+  a mapping of the space $X$ into itself. Given $[X, bb(Z)^+, f]$, $f(dot, 1)$
+  is its _transition function_. Thus if $[X, bb(Z)^+, f]$ is reversible, its
+  transition function is invertible, with inverse given by $f(dot, -1)$.
+
+  The set $f(p, S) = { f(p, i) | i in S }$ is called the _orbit_ or _trajectory_
+  of the point $p$. A _fixed point_ of $[X, S, f]$ is a point $p$ such that
+  $f(p, s) = p$ for all $s in S$. A set $A subset X$ is _invariant w.r.t._ $f$,
+  or simply _invariant_, if $f(A, s) subset A$ for all $s in S$.
+
+  The notions of equivalence and homomorphism are crucial. Two dynamical systems
+  $[X, S, f]$, $[Y, S, g]$ are said to be _isomorphic_ (also
+  _topologically equivalent_ or simply _equivalent_) if there exists a
+  homeomorphism $psi : X -> Y$ such that
+  $ psi(f(p, s)) = g(psi(p), s) $
+  for all $p in X$ and $s in S$. If the mapping $psi$ is only continuous, then
+  $[X, S, f]$ is said to be _homomorphic_ to $[Y, S, g]$. Homomorphisms preserve
+  trajectories, fixed points, and invariant sets.
+
+  In this paper, the continuous dynamical systems dealt with are defined by the
+  solutions of ordinary differential equations (ODEs):
+  $ x'(t) = f(x(t)) $
+  where $x(t) in X subset bb(R)^n$. The function $f : X -> bb(R)^n$ is called a
+  _vector field_ on $bb(R)^n$. The resulting dynamical system is then given by
+  $phi(x_0, t) = x(t)$ where $x(dot)$ is the solution to the above equation
+  starting at $x_0$ at $t = 0$. (We assume existence and uniqueness of
+  solutions; see [Hirsch, Smalle] for conditions). A system of ODEs is called
+  _autonomous_ or _time-invariant_ if its vector field does not depend
+  explicitly on time. Throughout, the shorthand _continuous_ (resp. _Lipschitz_)
+  _ODEs_ denotes ODEs with continuous (resp. Lipschitz) vector fields.
+
+  An ODE _with input and outputs_ is given by
+  $ x'(t) & = f(x(t), u(t)), \ y(t) & = h(x(t)) $
+  where $x(t) in X subset bb(R)^n$, $u(t) in U subset bb(R)^m$,
+  $y in Y subset bb(R)^p$, $f : bb(R)^n times bb(R)^m -> bb(R)^n$, and
+  $h : bb(R)^n -> bb(R)^p$. The functions $u(dot)$ and $y(dot)$ are the _inputs_
+  and _outputs_, respectively.
+])
+
 == Miscellaneous
 
 === Signals
@@ -507,7 +694,7 @@ $f : S u p e r d e n s e(bb(V))$, $f(t, n) = f(t + n partial)$ in the
 non-standard semantics.
 
 Our representation instead uses
-$S i g n a l(bb(V)) = S t r e a m((h : bb(R)) * ([0, h] -> bb(V)))$. Can we 
+$S i g n a l(bb(V)) = S t r e a m((h : bb(R)) * ([0, h] -> bb(V)))$. Can we
 convert between the two as we want?
 
 #breakl(```ocaml
diff --git a/doc/pres.typ b/doc/pres.typ
index 88c12a8..f4b52bd 100644
--- a/doc/pres.typ
+++ b/doc/pres.typ
@@ -145,7 +145,7 @@ Solution: simulate assertions with their own solver.
 Nodes take an additional reset parameter `'p` for their reset function:
 
 #align(top)[
-  #grid(columns: (3fr, 2fr), align: (left, left), 
+  #grid(columns: (3fr, 2fr), align: (left, left),
   ```ocaml
   type ('p, 'a, 'b) dnode =
     DNode : {
@@ -200,7 +200,7 @@ Hybrid nodes are a bit more complex:
 #slide(repeat: 3, self => [
   #let (uncover, only) = utils.methods(self)
 
-  ODE solvers are discrete nodes producing streams of functions defined on 
+  ODE solvers are discrete nodes producing streams of functions defined on
   successive intervals:
   $(h: bb(R)_+) -->^D (h': bb(R)_+) times (u: [0,h'] -> bb(V))$.
 
@@ -249,7 +249,7 @@ Hybrid nodes are a bit more complex:
 #slide(repeat: 3, self => [
   #let (uncover, only) = utils.methods(self)
 
-  Zero-crossing solvers are discrete nodes too, looking for zero-crossings on 
+  Zero-crossing solvers are discrete nodes too, looking for zero-crossings on
   functions:
   $(h: bb(R)_+) times (u: [0, h] -> bb(V)) -->^D
    (h': bb(R)_+) times (z: bb(Z)?)$.
@@ -456,7 +456,7 @@ When receiving a new value, store it
   box((10, 1.25), (3, 1.5), double: true, content: `state`)
 
   // (h, u) -> state
-  content((-2.25, 3), $(h, u)$) 
+  content((-2.25, 3), $(h, u)$)
   arrow((-1, 3), (r, 1.5), (d, 2.5), (r, 5), (u, 1.5), (r, 4.5))
 
   // sim -> bot
@@ -583,7 +583,7 @@ When receiving a new value, store it
 
   // solver.step
   box((6.5, 1), (3, 1.5), content: `step`)
-  
+
   // solver.state
   box((6.5, 2.75), (3, 1.5), double: true, content: `state`)
 
@@ -757,7 +757,7 @@ If no zero-crossing is found, there is no discontinuity.
     content((-3, 1.5), $"Signal"(alpha)$)
     content((5.5, 2), $"Signal"(beta)$)
     content((14, 1.5), $"Signal"(gamma)$)
-  
+
     content((-4, -1),   $alpha:$)
     content((-4, -2.5), $beta:$)
     content((-4, -4),   $gamma:$)
@@ -795,7 +795,7 @@ If no zero-crossing is found, there is no discontinuity.
   The output signal stops at least as often as the input signal does.
   #linebreak()
   This calls for a "lazy" composition: we request rather than provide data.
-  
+
   #align(center, cetz-canvas({
     import cetz.draw: *
     box((0, 0), (3, 3), content: $M_1$)
@@ -809,7 +809,7 @@ If no zero-crossing is found, there is no discontinuity.
     content((-3, 1.5), $"Signal"(alpha)$)
     content((5.5, 2), $"Signal"(beta)$)
     content((14, 1.5), $"Signal"(gamma)$)
-  
+
     content((-4, -1), $alpha:$)
     content((-4, -2.5), $beta:$)
     content((-4, -4), $gamma:$)
diff --git a/doc/rep.typ b/doc/rep.typ
new file mode 100644
index 0000000..fe4061f
--- /dev/null
+++ b/doc/rep.typ
@@ -0,0 +1,192 @@
+#import "@preview/cetz:0.4.0"
+#import "@preview/finite:0.4.1"
+
+#set page(paper: "a4")
+#set par(justify: true)
+#set raw(syntaxes: "zelus.sublime-syntax")
+#show raw: set text(font: "CaskaydiaCove NF")
+
+#let lustre = smallcaps[Lustre]
+#let esterel = smallcaps[Esterel]
+#let simulink = smallcaps[Simulink]
+#let sundials = smallcaps[Sundials CVODE]
+#let zelus = smallcaps[Zélus]
+#let TODO(..what) = {
+  let msg = if what.pos().len() == 0 { "" } else { ": " + what.pos().join("") }
+  block(fill: red, width: 100%, inset: 5pt, align(center, raw("TODO" + msg)))
+}
+#let adot(s) = $accent(#s, dot)$
+#let addot(s) = $accent(#s, dot.double)$
+
+#align(center)[
+  #text(16pt)[
+    *Internship Report --- MPRI M2 \
+    Hybrid System Models with Transparent Assertions*
+  ]
+
+  Henri Saudubray, supervised by Marc Pouzet, Inria PARKAS
+]\
+
+#heading(outlined: false)[General Context]
+What is the report about? Where does the problem come from? What is the state of
+the art in that area?
+
+Hybrid systems modelers such as #simulink are essential tools in the development
+of embedded systems evolving in physical environments. They allow for precise
+descriptions of hybrid systems through both continuous-time behaviour defined
+through Ordinary Differential Equations (ODEs) and discrete-time reactive
+behaviour similar to what is found in the synchronous languages such as #lustre.
+The #zelus language aims to reconcile synchronous languages with hybrid systems,
+by taking a synchronous language kernel and extending it with continuous-time
+constructs similar to what is found in tools like #simulink. Continuous-time
+behaviour is computed through the use of an off-the-shelf ODE solver such as
+#sundials.
+
+#heading(outlined: false)[Research Problem]
+What specific question(s) have you studied? What motivates the question? What
+are the applications/consequences? Is it a new problem? Why did you choose this
+problem?
+
+The simulation of hybrid system models, as done in #simulink and #zelus, uses a
+single ODE solver instance to simulate the entire model. This has various
+advantages: it is less computationally intensive, and simplifies the work of the
+compiler. Unfortunately, it also raises a difficult problem: sub-systems which
+seemingly should not interfere with each other end up affecting each other's
+results. This is due to the chosen integration method: an adaptive solver like
+#sundials will vary its step length throughout the integration process, and the
+addition of ODEs in the overall system can influence these step lengths,
+affecting the results obtained for pre-existing ODEs. This is particularly
+problematic in the case of runtime assertions, which are typically expected to
+be transparent: they should not affect the final result of the computation.
+
+#heading(outlined: false)[Proposed Contributions]
+What is your answer to the research problem? Please remain at a high level;
+technical details should be given in the main body of the report. Pay special
+attention to the description of the _scientific_ approach.
+
+To solve this, we propose a new runtime for the #zelus language, which simulates
+assertions with their own solvers in order to maintain the separation between
+assertions and the model they operate on.
+
+#TODO()
+
+#heading(outlined: false)[Arguments Supporting Their Validity]
+What is the evidence that your solution is a good solution? (Experiments?
+Proofs?) Comment on the robustness of your solution: does it rely on
+assumptions, and if so, to what extent?
+
+#TODO("Justify")
+
+#heading(outlined: false)[Summary and Future Work]
+What did you contribute to the area? What comes next? What is a good _next_ step
+or question?
+
+#TODO("Next steps")
+
+#pagebreak()
+#outline()
+#pagebreak()
+
+= Background
+
+== Hybrid systems
+
+Hybrid systems are dynamical systems whose behaviour evolves through both
+continuous and discrete phases. They are used to model software interacting with
+physical systems. Continuous phases are described using ordinary differential
+equations (ODEs), while discrete phases can be represented as a reactive
+program in a synchronous language such as #lustre or #esterel.
+
+As a first example, say we wish to model a bouncing ball. We could start by
+describing its distance from the ground $y$ with a second-order differential
+equation
+$ addot(y) = -9.81 $
+where $addot(y)$ denotes the second order derivative of $y$ with
+respect to time (the acceleration of the ball), and $9.81$ is the gravitational
+constant $g$: the acceleration of the ball is its negation. We now give the
+initial position and speed of the ball:
+$ y(0) = 50 space space space adot(y)(0) = 0 $
+We have just described an initial value problem: given an ODE and an initial
+value for its dependent variable, its solution is a function $y(t)$ returning
+the value of the variable at a given time $t$. We can find an approximation of
+this function using an ODE solver, such as #sundials.
+
+As of right now, our ball will fall until the end of time; we have not said
+anything about how it bounces when it hits the floor. To do so, we need a notion
+of discrete _events_. These are modelled by zero-crossings: we monitor a certain
+value and stop when it goes from strictly negative to positive or null. For our
+purposes, we choose $-y$ as the monitored value, and call the zero-crossing
+event $z$. When $z$ occurs (i.e., when the ball touches the ground), we set the
+speed $adot(y)$ to $-k dot "last"(adot(y))$, where $"last"(y)$ denotes the left
+limit of $y$ (we cannot specify $adot(y)$ in terms of itself), and $k$ is a
+factor modelling the loss of inertia due to the collision (say, $0.8$). We can
+then resume the approximation of the solution.
+
+@lst:ball.zls shows how such a model might be expressed in the concrete syntax
+of #zelus.
+
+#figure(placement: top, gap: 2em,
+  ```zelus
+  let hybrid ball () = y where
+    rec der y = y' init 50.0
+    and der y' = -9.81 init 0.0 reset z -> -0.8 *. (last y')
+    and z = up (-. y)
+  ```,
+  caption: [The bouncing ball in #zelus]
+) <lst:ball.zls>
+
+More formally, a hybrid system can be described as an automaton
+
+== Executing models
+
+Executing such a model is quite simple. There are two modes of execution:
+discrete and continuous. In continuous mode, we call the ODE solver in order to
+obtain an approximation of the variables defined through ODEs, and monitor for
+zero-crossings. If a zero-crossing occurs, we enter the discrete mode, in which
+we perform computation steps as needed, until no other zero-crossing occurs, in
+which case we go back to the continuous mode, and repeat, as seen in @automaton.
+
+#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)
+), caption: [High-level overview of the runtime], placement: top) <automaton>
+
+= Runtime
+To solve this issue, we need to redefine what the runtime of our hybrid system
+models is. The runtime is the piece of software that handles the pairing of a
+hybrid model with a solver and the different execution phases that need to take
+place.
+
+To allow for independent simulation of the assertions, we need to simulate them
+with their own ODE solvers. Since assertions can themselves contain ODEs, they
+can be viewed as hybrid systems themselves. Assertions can also contain their
+own sub-assertions, and so recursively: our runtime needs to handle this
+accordingly.
+
+== Hybrid models with assertions
+A hybrid model with assertions can be defined using a recursive data type:
+```ocaml
+type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a =
+  HNodeA : {
+    body : ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec;
+    assertions : ('p, 's, unit, 'y, 'yder, 'zin, 'zout) hnode_a list;
+  } -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a
+```
+
+Notice that the list of assertions takes as input the inner state of the parent
+model: assertions are checked on the model state, and any variable which is
+required by the assertion becomes a state variable.
+
+== Solvers as synchronous nodes
+== Simulations as synchronous nodes
+#TODO("talk about the new runtime")
+
+= Assertions
+#TODO("talk about how assertions are done")
+
+#pagebreak()
+= Annex
+#set heading(level: 2)
+#bibliography("sources.bib", full: true)
+#set heading(level: 1)
diff --git a/exm/ball.ml b/exm/ball.ml
index dcdd63c..d78109d 100644
--- a/exm/ball.ml
+++ b/exm/ball.ml
@@ -29,7 +29,7 @@ let fder y yd =
   else begin yd.{0} <- 0.0; yd.{1} <- 0.0; yd.{2} <- 0.0; yd.{3} <- 0.0 end;
   yd
 let fzer y zo  = zo.{0} <- -. y.{1}; zo
-let fout _ _ y = of_array [| y.{1} |]
+let fout _ _ y = of_array [| y.{1}; y.{0} |]
 let jump _     = true
 let horizon _  = max_float
 let cget s     = s.lx
@@ -38,9 +38,9 @@ let zset s zin = { s with zin }
 let step ({ zin; lx; _ } as s) zfalse =
   let lx = if zin.{0} = 1l then
     of_array [| -. 0.8 *. lx.{0}; 0.0; lx.{2}; lx.{3} |] else lx in
-  of_array [| s.lx.{1} |], { zin=zfalse; lx; i=false }
+  of_array [| s.lx.{1}; s.lx.{0} |], { zin=zfalse; lx; i=false }
 
-let bouncing_ball () : (_, _, carray, carray, carray, zarray, carray) hnode =
+let bouncing_ball () : (_, _, _, carray, carray, carray, zarray, carray) hrec =
   let yd         = cmake csize in
   let zout       = cmake zsize in
   let zfalse     = zmake 1 in
@@ -49,10 +49,10 @@ let bouncing_ball () : (_, _, carray, carray, carray, zarray, carray) hnode =
   let step s _   = step s zfalse in
   let state      = { zin=zfalse; lx=of_array [|y'0;y0;x'0;x0|]; i=true } in
   let reset _ _  = state in
-  HNode { state; fder; fzer; fout; step; reset; horizon;
-          jump; cset; cget; zset; csize; zsize }
+  { state; fder; fzer; fout; step; reset; horizon; jump; cset; cget; zset;
+    csize; zsize }
 
 let errmsg = "Too many arguments for the model (needed: 0)"
 let init = function
-  | [] -> bouncing_ball ()
+  | [] -> HNode (bouncing_ball ())
   | _  -> raise (Invalid_argument errmsg)
diff --git a/exm/dune b/exm/dune
index 9eb098f..aa40c69 100644
--- a/exm/dune
+++ b/exm/dune
@@ -1,3 +1,10 @@
+(env
+ (dev
+  (flags
+   (:standard -w -a))))
+
 (library
-  (name examples)
-  (libraries hsim solvers))
+ (name examples)
+ (libraries hsim solvers))
+
+(include_subdirs unqualified)
diff --git a/exm/sin1x.ml b/exm/sin1x.ml
new file mode 100644
index 0000000..a6f2447
--- /dev/null
+++ b/exm/sin1x.ml
@@ -0,0 +1,72 @@
+
+(* Example due to JB. Jeannin - zero-crossing detection *)
+(* the signal x(t) is expected to have a zero-crossing at [t = 0.3] *)
+
+(* to draw: ./main.exe -s sin_1_x -stop 2 -sample 10000 | feedgnuplot --stream --domain --lines *)
+
+open Hsim.Types
+open Solvers.Zls
+
+(* let hybrid sin_1_x() =
+     let der t = 1.0 init 0.0 in
+     let t = t -. 0.3 in
+     let v = sin(1/t) in
+     let o = t *. (v *. v) in [that is: t *. (sin(1/t)^2]
+     present up(o) -> 1.0 else 0.0, o *)
+
+let of_array a = Bigarray.Array1.of_array Bigarray.Float64 Bigarray.c_layout a
+
+type ('a, 'b) state = { s_x: 'a; zin: 'b }
+
+let csize = 1
+let zsize = 1
+
+let fder _ _ yd = yd.{0} <- 1.0
+
+let fzer _ y zout =
+  let t = y.{0} in
+  let t = t -. 0.3 in
+  let v = sin(1.0 /. t) in
+  let o = t *. (v *. v) in
+  zout.{0} <- o
+
+let fout s y =
+  let z = if s.zin.{0} = 1l then 1.0 else 0.0 in
+  let t = y.{0} in
+  let t = t -. 0.3 in
+  let v = sin(1.0 /. t) in
+  let o = t *. (v *. v) in
+  of_array [| z; o |]
+
+let step ({ s_x; zin } as s) zfalse =
+  let z = if zin.{0} = 1l then 1.0 else 0.0 in
+  let t = s_x.{0} in
+  let t = t -. 0.3 in
+  let v = sin(1.0 /. t) in
+  let o = t *. (v *. v) in
+  of_array [| z; o |], { s with zin = zfalse }
+
+let reset _ s = s
+
+let cget { s_x; _ } = s_x
+let cset s l_x = { s with s_x = l_x }
+let zset s zin = { s with zin }
+let jump _ = true
+let horizon _ = max_float
+
+let sin_1_x () =
+  let zfalse = zmake 1 in
+  let yd     = cmake 1 in
+  let zout   = cmake 1 in
+  let fder _ _ y = fder 0.0 y yd; yd in
+  let fzer _ _ y = fzer 0.0 y zout; zout in
+  let fout s _ y = fout s y in
+  let step s _ = step s zfalse in
+  let state = { s_x = of_array [| 0.0 |]; zin = zfalse } in
+  { state; fder; fzer; fout; step; reset; horizon;
+    cset; cget; zset; csize; zsize; jump }
+
+let errmsg = "Too many arguments for the model (needed: 0)"
+let init = function
+  | [] -> HNode (sin_1_x ())
+  | _  -> raise (Invalid_argument errmsg)
diff --git a/exm/sin1x_der.ml b/exm/sin1x_der.ml
new file mode 100644
index 0000000..e4bf116
--- /dev/null
+++ b/exm/sin1x_der.ml
@@ -0,0 +1,79 @@
+
+(* Example due to JB. Jeannin - zero-crossing detection. *)
+
+(* The same as [sin1x.ml] but with an integrator.
+
+   d(t . (sin(1/t))^2))/dt = sin(1/t)(sin(1/t) - (2 . cos(1/t)) / t)
+
+   d(sin(f(x)))/dx = cos(f(x)) f'(x) *)
+
+open Hsim.Types
+open Solvers.Zls
+
+(* let hybrid sin_1_x() =
+     let der t0 = 1.0 init 0.0 in
+     let d = 0.6 in
+     let t = t0 -. d in
+     let der o = sin(1.0 /. t) *. (sin(1.0 /. t) -. (2.0 *. cos(1.0 /. t)) /. t)
+         init (-. d *. (sin(1.0 /. (-. d))) ** 2.0) in
+     present up(o) -> 1.0 else 0.0, o *)
+
+let of_array a = Bigarray.Array1.of_array Bigarray.Float64 Bigarray.c_layout a
+
+type ('a, 'b) state = { s_x: 'a; zin: 'b }
+
+let d = 0.66
+let y0 = -. d *. let v = (sin(1.0 /. (-. d))) in v *. v
+let csize = 2
+let zsize = 1
+
+let fder d y yd =
+  yd.{0} <- 1.0;
+  let t = y.{0} -. d in
+  yd.{1} <- sin(1.0 /. t) *. (sin(1.0 /. t) -. (2.0 *. cos(1.0 /. t)) /. t)
+
+let fzer z y zout =
+  let o = y.{1} in
+  zout.{0} <- if z then o else 1.0
+
+let fout s y =
+  let z = if s.zin.{0} = 1l then 1.0 else 0.0 in
+  let o = y.{1} in
+  of_array [| z; o |]
+
+let step ({ s_x; zin } as s) zfalse =
+  let z = if zin.{0} = 1l then 1.0 else -1.0 in
+  let o = s_x.{1} in
+  of_array [| z; o |], { s with zin = zfalse }
+
+let reset _ s = s
+
+let cget { s_x; _ } = s_x
+let cset s l_x = { s with s_x = l_x }
+let zset s zin = { s with zin }
+let horizon _ = max_float
+let jump _ = true
+
+let sin_1_x z d =
+  let zfalse = zmake 1 in
+  let yd     = cmake 2 in
+  let zout   = cmake 1 in
+  let fder _ _ y = fder d y yd; yd in
+  let fzer _ _ y = fzer z y zout; zout in
+  let fout s _ y = fout s y in
+  let step s _ = step s zfalse in
+  let state = { s_x = of_array [| 0.0; y0 |]; zin = zfalse } in
+  { state; fder; fzer; fout; step; horizon;
+    cset; cget; zset; csize; zsize; reset; jump }
+
+let errmsg_many    = "Too many arguments to model (needed: [bool, float])"
+let errmsg_few     = "Too few arguments to model (needed: [bool, float])"
+let errmsg_invalid = "Invalid arguments to model (needed: [bool, float])"
+
+let init = function
+  | [zer; d] ->
+      let zer, d = try bool_of_string zer, float_of_string d
+      with _ -> raise (Invalid_argument errmsg_invalid) in
+      HNode (sin_1_x zer d)
+  | [] | [_] -> raise (Invalid_argument errmsg_few)
+  | _  -> raise (Invalid_argument errmsg_many)
diff --git a/exm/sincos.ml b/exm/sincos.ml
index 17131c6..ab62149 100644
--- a/exm/sincos.ml
+++ b/exm/sincos.ml
@@ -34,11 +34,11 @@ let sinus_cosinus theta0 omega =
   HNode { state; fder; fzer; fout; step; reset; horizon;
           jump; cset; cget; zset; csize; zsize }
 
-let errmsg_invalid = "Invalid arguments to model (needed: 2 floats)"
-let errmsg_few     = "Too few arguments to model (needed: 2 floats)"
-let errmsg_many    = "Too many arguments to model (needed: 2 floats)"
+let errmsg_invalid = "Invalid arguments to model (needed: [float, float])"
+let errmsg_few     = "Too few arguments to model (needed: [float, float])"
+let errmsg_many    = "Too many arguments to model (needed: [float, float])"
 let init = function
-  | [t0; om] ->
+  | [om; t0] ->
       let t0, om = try float_of_string t0, float_of_string om
       with Failure _ -> raise (Invalid_argument errmsg_invalid) in
       sinus_cosinus t0 om
diff --git a/exm/zelus/ballz/ballz.zls b/exm/zelus/ballz/ballz.zls
new file mode 100644
index 0000000..5f7e95f
--- /dev/null
+++ b/exm/zelus/ballz/ballz.zls
@@ -0,0 +1,22 @@
+let g = 9.81
+let y0 = 50.0
+let y'0 = 0.0
+
+let hybrid ball (y0, y'0) = (y, y', z) where
+  rec der y  = y'   init y0
+  and der y' = -. g init y'0 reset z -> -0.8 *. (last y')
+  and z = up(-. y)
+
+let hybrid main () =
+  let der t = 1.0 init 0.0 in
+  let rec der p = 1.0 init -0.01 reset s -> -0.01
+      and s = up(p) in 
+  let (y, y', z) = ball (y0, y'0) in
+  present z | s -> (
+    print_float t;
+    print_string "\t";
+    print_float y;
+    print_string "\t";
+    print_float y';
+    print_newline ()
+  ); ()
diff --git a/exm/zelus/ballz/dune b/exm/zelus/ballz/dune
new file mode 100644
index 0000000..81087c1
--- /dev/null
+++ b/exm/zelus/ballz/dune
@@ -0,0 +1,6 @@
+(rule
+ (targets ballz.ml ballz.zci)
+ (deps
+  (:zl ballz.zls))
+ (action
+  (run zeluc %{zl})))
diff --git a/exm/zelus/sin_1_x.ml b/exm/zelus/sin_1_x.ml
new file mode 100644
index 0000000..d185ae5
--- /dev/null
+++ b/exm/zelus/sin_1_x.ml
@@ -0,0 +1,74 @@
+(* The Zelus compiler, version 2.2-dev
+  (2025-06-16-15:24) *)
+open Common
+open Ztypes
+open Solvers
+(* open Zls *)
+
+type ('f , 'e , 'd , 'c , 'b , 'a) _sin_1_x =
+  { mutable major_22 : 'f ;
+    mutable i_29 : 'e ;
+    mutable x_28 : 'd ;
+    mutable result_27 : 'c ; mutable o_26 : 'b ; mutable t0_23 : 'a }
+
+let sin_1_x (cstate_30:Ztypes.cstate) =
+
+  let sin_1_x_alloc _ =
+    cstate_30.cmax <- (+) cstate_30.cmax  2 ;
+    cstate_30.zmax <- (+) cstate_30.zmax  1;
+    { major_22 = false ;
+      i_29 = (false:bool) ;
+      x_28 = { zin = false; zout = 1. } ;
+      result_27 = (42.:float) ;
+      o_26 = { pos = 42.; der = 0. } ; t0_23 = { pos = 42.; der = 0. } } in
+  let sin_1_x_step self ((_time_21:float) , ()) =
+    ((let (cindex_31:int) = cstate_30.cindex in
+      let cpos_33 = ref (cindex_31:int) in
+      let (zindex_32:int) = cstate_30.zindex in
+      let zpos_34 = ref (zindex_32:int) in
+      cstate_30.cindex <- (+) cstate_30.cindex  2 ;
+      cstate_30.zindex <- (+) cstate_30.zindex  1 ;
+      self.major_22 <- cstate_30.major ;
+      (if cstate_30.major then
+       for i_1 = cindex_31 to 1 do Zls.set cstate_30.dvec  i_1  0. done
+       else ((self.o_26.pos <- Zls.get cstate_30.cvec  !cpos_33 ;
+              cpos_33 := (+) !cpos_33  1) ;
+             (self.t0_23.pos <- Zls.get cstate_30.cvec  !cpos_33 ;
+              cpos_33 := (+) !cpos_33  1))) ;
+      (let (result_35:(float  * float)) =
+           (if self.i_29 then
+            self.o_26.pos <- ( *. ) ((~-.) 0.6)
+                                    (( ** ) (sin ((/.) 1.  ((~-.) 0.6)))  2.))
+           ;
+           self.i_29 <- false ;
+           self.t0_23.der <- 1. ;
+           (let (t_25:float) = (-.) self.t0_23.pos  0.6 in
+            self.o_26.der <- ( *. ) (sin ((/.) 1.  t_25))
+                                    ((-.) (sin ((/.) 1.  t_25))
+                                          ((/.) (( *. ) 2.
+                                                        (cos ((/.) 1.  t_25)))
+                                                 t_25)) ;
+            self.x_28.zout <- self.o_26.pos ;
+            (begin match self.x_28.zin with
+                   | true -> self.result_27 <- 1.
+                   | _ -> self.result_27 <- 0.  end) ;
+            (self.result_27 , self.o_26.pos)) in
+       cpos_33 := cindex_31 ;
+       (if cstate_30.major then
+        (((Zls.set cstate_30.cvec  !cpos_33  self.o_26.pos ;
+           cpos_33 := (+) !cpos_33  1) ;
+          (Zls.set cstate_30.cvec  !cpos_33  self.t0_23.pos ;
+           cpos_33 := (+) !cpos_33  1)) ; ((self.x_28.zin <- false)))
+        else (((self.x_28.zin <- Zls.get_zin cstate_30.zinvec  !zpos_34 ;
+                zpos_34 := (+) !zpos_34  1)) ;
+              zpos_34 := zindex_32 ;
+              ((Zls.set cstate_30.zoutvec  !zpos_34  self.x_28.zout ;
+                zpos_34 := (+) !zpos_34  1)) ;
+              ((Zls.set cstate_30.dvec  !cpos_33  self.o_26.der ;
+                cpos_33 := (+) !cpos_33  1) ;
+               (Zls.set cstate_30.dvec  !cpos_33  self.t0_23.der ;
+                cpos_33 := (+) !cpos_33  1)))) ; result_35)):float * float) in
+
+  let sin_1_x_reset self  =
+    ((self.i_29 <- true ; self.t0_23.pos <- 0.):unit) in
+  Node { alloc = sin_1_x_alloc; step = sin_1_x_step ; reset = sin_1_x_reset }
diff --git a/exm/zelus/sin_1_x.zci b/exm/zelus/sin_1_x.zci
new file mode 100644
index 0000000..42b4e83
Binary files /dev/null and b/exm/zelus/sin_1_x.zci differ
diff --git a/exm/zelus/sin_1_x.zls b/exm/zelus/sin_1_x.zls
new file mode 100644
index 0000000..48cc71e
--- /dev/null
+++ b/exm/zelus/sin_1_x.zls
@@ -0,0 +1,7 @@
+let hybrid sin_1_x () =
+  let der t0 = 1.0 init 0.0 in
+  let d = 0.6 in
+  let t = t0 -. d in
+  let der o = sin(1.0 /. t) *. (sin(1.0 /. t) -. (2.0 *. cos(1.0 /. t)) /. t)
+			  init (-. d *. (sin(1.0 /. (-. d))) ** 2.0) in
+  (present up(o) -> 1.0 else 0.0), o
diff --git a/exm/zelus/sincos/dune b/exm/zelus/sincos/dune
new file mode 100644
index 0000000..ef360bb
--- /dev/null
+++ b/exm/zelus/sincos/dune
@@ -0,0 +1,6 @@
+(rule
+ (targets sincosz.ml sincosz.zci)
+ (deps
+  (:zl sincosz.zls))
+ (action
+  (run zeluc %{zl})))
diff --git a/exm/zelus/sincos/sincosz.zls b/exm/zelus/sincos/sincosz.zls
new file mode 100644
index 0000000..bfaa71d
--- /dev/null
+++ b/exm/zelus/sincos/sincosz.zls
@@ -0,0 +1,4 @@
+
+let hybrid f () = (sin, cos) where
+  rec der sin = cos init 0.0
+  and der cos = -. sin init 1.0
diff --git a/exm/ztypes.ml b/exm/ztypes.ml
new file mode 100644
index 0000000..567f8ff
--- /dev/null
+++ b/exm/ztypes.ml
@@ -0,0 +1,4 @@
+include Common
+include Ztypes
+include Solvers
+
diff --git a/src/bin/lift.ml b/src/bin/lift.ml
new file mode 100644
index 0000000..ca9f336
--- /dev/null
+++ b/src/bin/lift.ml
@@ -0,0 +1,118 @@
+
+open Hsim.Types
+open Solvers.Zls
+open Common.Ztypes
+
+type ('s, 'a) state =
+  { mutable state : 's; mutable input : 'a option; mutable time : time }
+
+let lift f =
+  let cstate =
+    { cvec = cmake 0; dvec = cmake 0; cindex = 0; zindex = 0;
+      cend = 0; zend = 0; cmax = 0; zmax = 0;
+      zinvec = zmake 0; zoutvec = cmake 0;
+      major = false; horizon = max_float } in
+  let Node { alloc=f_alloc; step=f_step; reset=f_reset } = f cstate in
+  let state = { state = f_alloc (); input = None; time = 0.0 } in
+  let csize, zsize = cstate.cmax, cstate.zmax in
+  let no_roots_in  = zmake zsize in
+  let no_roots_out = cmake zsize in
+  let ignore_der   = cmake csize in
+  let cstates      = cmake csize in
+  cstate.cvec <- cstates;
+  f_reset state.state;
+
+  let no_time = -1.0 in
+
+  (* the function that compute the derivatives *)
+  let fder { state; _ } input y =
+    cstate.major <- false;
+    cstate.zinvec <- no_roots_in;
+    cstate.zoutvec <- no_roots_out;
+    cstate.cvec <- y;
+    cstate.dvec <- ignore_der;
+    cstate.cindex <- 0;
+    cstate.zindex <- 0;
+    ignore (f_step state (no_time, input));
+    cstate.dvec in
+
+  (* the function that compute the zero-crossings *)
+  let fzer { state; _ } input y =
+    cstate.major <- false;
+    cstate.zinvec <- no_roots_in;
+    cstate.dvec <- ignore_der;
+    cstate.zoutvec <- no_roots_out;
+    cstate.cvec <- y;
+    cstate.cindex <- 0;
+    cstate.zindex <- 0;
+    ignore (f_step state (no_time, input));
+    cstate.zoutvec in
+
+  (* the function which compute the output during integration *)
+  let fout { state; _ } input y =
+    cstate.major <- false;
+    cstate.zoutvec <- no_roots_out;
+    cstate.dvec <- ignore_der;
+    cstate.zinvec <- no_roots_in;
+    cstate.cvec <- y;
+    cstate.cindex <- 0;
+    cstate.zindex <- 0;
+    f_step state (no_time, input) in
+
+  (* the function which compute a discrete step *)
+  let step ({ state; time; _ } as st) input =
+    cstate.major <- true;
+    cstate.horizon <- infinity;
+    cstate.zinvec <- no_roots_in;
+    cstate.zoutvec <- no_roots_out;
+    cstate.dvec <- ignore_der;
+    cstate.cindex <- 0;
+    cstate.zindex <- 0;
+    let o = f_step state (time, input) in
+    o, { st with state; input=Some input } in
+
+  let reset _ ({ state; _ } as st) = f_reset state; st in
+
+  (* horizon *)
+  let horizon _ = cstate.horizon in
+
+  let jump _ = true in
+
+  (* the function which sets the zinvector into the *)
+  (* internal zero-crossing variables *)
+  let zset ({ state; input; _ } as st) zinvec =
+    cstate.major <- false;
+    cstate.zoutvec <- no_roots_out;
+    cstate.dvec <- ignore_der;
+    cstate.zinvec <- zinvec;
+    cstate.cindex <- 0;
+    cstate.zindex <- 0;
+    ignore (f_step state (no_time, Option.get input));
+    st in
+
+  let cset ({ state; input; _ } as st) _ =
+    cstate.major <- false;
+    cstate.horizon <- infinity;
+    cstate.zinvec <- no_roots_in;
+    cstate.zoutvec <- no_roots_out;
+    cstate.dvec <- ignore_der;
+    cstate.cindex <- 0;
+    cstate.zindex <- 0;
+    ignore (f_step state (no_time, Option.get input));
+    st in
+
+  let cget { state; input; _ } =
+    cstate.major <- false;
+    cstate.horizon <- infinity;
+    cstate.zinvec <- no_roots_in;
+    cstate.zoutvec <- no_roots_out;
+    cstate.dvec <- ignore_der;
+    cstate.cindex <- 0;
+    cstate.zindex <- 0;
+    ignore (f_step state (no_time, Option.get input));
+    cstate.cvec in
+
+  HNode
+    { state; fder; fzer; step; fout; reset;
+      horizon; cset; cget; zset; zsize; csize; jump }
+
diff --git a/src/bin/main.ml b/src/bin/main.ml
index 983c89d..3d6c608 100644
--- a/src/bin/main.ml
+++ b/src/bin/main.ml
@@ -4,17 +4,26 @@ open Solvers
 open Examples
 open Common
 open Types
+open Lift
 
-let sample   = ref 10
-let stop     = ref 30.0
-let greedy   = ref false
+let sample   = ref 1
+let stop     = ref 10.0
+let accel    = ref false
 let inplace  = ref false
 let sundials = ref false
+let speed    = ref false
 let steps    = ref 1
 let model    = ref None
+let minstep  = ref None
+let maxstep  = ref None
+let mintol   = ref None
+let maxtol   = ref None
+let no_print = ref false
+let zelus    = ref false
 
 let gt0i v i  = v := if i <= 0   then 1   else i
 let gt0f v f  = v := if f <= 0.0 then 1.0 else f
+let opt r s = r := Some s
 let modelargs = ref []
 
 let set_model s =
@@ -26,46 +35,78 @@ 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 (doesn't support -accelerate)";
   "-inplace", Arg.Set inplace, "\tUse imperative solvers";
   "-steps", Arg.Int (gt0i steps), "n \tSplit into [n] steps (default=1)";
+  "-speed", Arg.Set speed, "\tLog the step length";
+  "-minstep", Arg.String (opt minstep), "\tSet minimum solver step length";
+  "-maxstep", Arg.String (opt maxstep), "\tSet maximum solver step length";
+  "-mintol",  Arg.String (opt mintol),  "\tSet minimum solver tolerance";
+  "-maxtol",  Arg.String (opt maxtol),  "\tSet maximum solver tolerance";
+  "-no-print", Arg.Set no_print, "\tDo not print output values";
+  "-zelus", Arg.Set zelus, "\tUse the output of the Zélus compiler";
 ]
 
 let errmsg = "Usage: " ^ Sys.executable_name ^ " [OPTIONS] MODEL\nOptions are:"
 
 let () = try Arg.parse (Arg.align opts) set_model errmsg with _ -> exit 2
 
-let m =
-  try match !model with
-  | None          -> Format.eprintf "Missing model\n"; exit 2
-  | Some "ball"   -> Ball.init !modelargs
-  | Some "vdp"    -> Vdp.init !modelargs
-  | Some "sincos" -> Sincos.init !modelargs
-  | Some "sqrt"   -> Sqrt.init !modelargs
-  | Some s        -> Format.eprintf "Unknown model: %s\n" s; exit 2
-  with Invalid_argument s -> Format.eprintf "%s\n" s; exit 2
+let args = List.rev !modelargs
+let () = ignore lift
 
-let z = StatefulZ.(if !inplace then InPlace.zsolve else Functional.zsolve)
+let wrap_zelus (HNode m) =
+  let ret = Bigarray.(Array1.create Float64 c_layout 0) in
+  let fout s a y = ignore (m.fout s a y); ret in
+  let step s () = let _, s = m.step s () in ret, s in
+  HNode { m with fout; step }
+
+let m =
+  try
+    if !zelus then
+      match !model with
+      | None           -> Format.eprintf "Missing model\n"; exit 2
+      | Some "ballz"   -> wrap_zelus (lift Ballz.main)
+      | Some "sincosz" -> wrap_zelus (lift Sincosz.f)
+      | Some s         -> Format.eprintf "Unknown model: %s\n" s; exit 2
+    else
+      match !model with
+      | None          -> Format.eprintf "Missing model\n"; exit 2
+      | Some "ball"   -> Ball.init args
+      | Some "vdp"    -> Vdp.init args
+      | Some "sincos" -> Sincos.init args
+      | Some "sqrt"   -> Sqrt.init args
+      | Some "sin1x"  -> Sin1x.init args
+      | Some "sin1xd" -> Sin1x_der.init args
+      | Some s        -> Format.eprintf "Unknown model: %s\n" s; exit 2
+  with Invalid_argument s -> Format.eprintf "%s\n" s; exit 2
 
 let st = if !inplace then (module State.InPlaceSimState    : State.SimState)
                      else (module State.FunctionalSimState : State.SimState)
 
+let output =
+  if !no_print then Output.ignore
+  else if !speed then Output.print_h
+  else Output.print (* Output.ignore *)
+
 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 open StatefulZ in
+    let z = if !inplace then InPlace.zsolve else Functional.zsolve in
+    let s = Solver.solver c (d_of_dc z) in
+    let open Sim.Sim(val st) in
+    run_until_n (output !sample (run m s))
   else
     let open StatefulRK45 in
     let c = if !inplace then InPlace.csolve else Functional.csolve in
+    let open StatefulZ in
+    let z = if !inplace then InPlace.zsolve else Functional.zsolve 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 !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 715a739..e60ea97 100644
--- a/src/bin/output.ml
+++ b/src/bin/output.ml
@@ -1,8 +1,8 @@
 
 open Hsim.Types
-open Common
+open Hsim.Utils
 
-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 +12,53 @@ let print_entry t y =
   Format.printf "\n";
   flush stdout
 
-let print samples { start; length; u } =
-  let step = length /. (float_of_int samples) in
+let print_entry_h y t h =
+  let n = Bigarray.Array1.dim y in
+  let rec loop i =
+    if i = n then ()
+    else (Format.printf "\t% .10e" y.{i}; loop (i+1)) in
+  Format.printf "% .10e\t% .10e" t h;
+  loop 0;
+  Format.printf "\n";
+  flush stdout
+
+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 (start +. length) (u length)
+    else if i = samples then print_entry (u h) (now +. h)
     else let t = float_of_int i *. step in
-    (print_entry (start +. t) (u t); loop (i+1)) in
-  if length <= 0.0 then begin Debug.print "D: "; print_entry start (u 0.0) end
-  else begin Debug.print "C: "; loop 0 end
+    (print_entry (u t) (now +. t); loop (i+1)) in
+  if h <= 0.0 then print_entry (u 0.0) now else loop 0
 
-let print_limits { start; length; _ } =
-  if length <= 0.0 then Format.printf "D: % .10e\n" start
-  else Format.printf "C: % .10e\t% .10e\n" start (start +. length)
+let print_sample_h 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) (now +. h) h
+    else let t = float_of_int i *. step in
+    (print_entry_h (u t) (now +. t) h; loop (i+1)) in
+  if h <= 0.0 then print_entry_h (u 0.0) now h else loop 0
+
+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, ((), ())) }
+
+let print_h samples n =
+  let DNode m = compose n (compose track (map (print_sample_h samples))) in
+  DNode { m with reset=fun p -> m.reset (p, ((), ())) }
+
+let ignore _ n =
+  let state = () in
+  let step () = function
+    | None -> None, ()
+    | Some _ -> Some (), () in
+  let reset () () = () in
+  let i = DNode { state; step; reset } in
+  let DNode n = compose n i in
+  DNode { n with reset=fun p -> n.reset (p, ()) }
 
diff --git a/src/lib/common/debug.ml b/src/lib/common/debug.ml
index e34b3cf..8db8349 100644
--- a/src/lib/common/debug.ml
+++ b/src/lib/common/debug.ml
@@ -2,3 +2,10 @@
 let debug = ref false
 
 let print s = if !debug then Format.printf "%s\n" s else ()
+
+let print_entry y =
+  let n = Bigarray.Array1.dim y in
+  let rec loop i =
+    if i = n then ()
+    else (Format.printf "\t% .10e" y.{i}; loop (i+1)) in
+  if !debug then (loop 0; Format.printf "\n"; flush stdout)
diff --git a/src/lib/common/ztypes.ml b/src/lib/common/ztypes.ml
new file mode 100644
index 0000000..82ee092
--- /dev/null
+++ b/src/lib/common/ztypes.ml
@@ -0,0 +1,151 @@
+(**************************************************************************)
+(*                                                                        *)
+(*                                Zelus                                   *)
+(*               A synchronous language for hybrid systems                *)
+(*                       http://zelus.di.ens.fr                           *)
+(*                                                                        *)
+(*                    Marc Pouzet and Timothy Bourke                      *)
+(*                                                                        *)
+(*  Copyright 2012 - 2019. All rights reserved.                           *)
+(*                                                                        *)
+(*  This file is distributed under the terms of the CeCILL-C licence      *)
+(*                                                                        *)
+(*  Zelus is developed in the INRIA PARKAS team.                          *)
+(*                                                                        *)
+(**************************************************************************)
+
+(* Type declarations and values that must be linked with *)
+(* the generated code *)
+type 'a continuous = { mutable pos: 'a; mutable der: 'a }
+
+type ('a, 'b) zerocrossing = { mutable zin: 'a; mutable zout: 'b }
+
+type 'a signal = 'a * bool
+type zero = bool
+
+(* a synchronous stream function with type 'a -D-> 'b *)
+(* is represented by an OCaml value of type ('a, 'b) node *)
+type ('a, 'b) node =
+    Node:
+      { alloc : unit -> 's; (* allocate the state *)
+        step : 's -> 'a -> 'b; (* compute a step *)
+        reset : 's -> unit; (* reset/inialize the state *)
+      } -> ('a, 'b) node
+
+(* the same with a method copy *)
+type ('a, 'b) cnode =
+    Cnode:
+      { alloc : unit -> 's; (* allocate the state *)
+        copy : 's -> 's -> unit; (* copy the source into the destination *)
+        step : 's -> 'a -> 'b; (* compute a step *)
+        reset : 's -> unit; (* reset/inialize the state *)
+      } -> ('a, 'b) cnode
+
+open Bigarray
+
+type time = float
+type cvec = (float, float64_elt, c_layout) Array1.t
+type dvec = (float, float64_elt, c_layout) Array1.t
+type zinvec = (int32, int32_elt,   c_layout) Array1.t
+type zoutvec = (float, float64_elt, c_layout) Array1.t
+
+(* The interface with the ODE solver *)
+type cstate =
+  { mutable dvec : dvec; (* the vector of derivatives *)
+    mutable cvec : cvec; (* the vector of positions *)
+    mutable zinvec : zinvec; (* the vector of boolean; true when the
+                                solver has detected a zero-crossing *)
+    mutable zoutvec : zoutvec; (* the corresponding vector that define
+                                  zero-crossings *)
+    mutable cindex : int; (* the position in the vector of positions *)
+    mutable zindex : int; (* the position in the vector of zero-crossings *)
+    mutable cend : int; (* the end of the vector of positions *)
+    mutable zend : int; (* the end of the zero-crossing vector *)
+    mutable cmax : int; (* the maximum size of the vector of positions *)
+    mutable zmax : int; (* the maximum number of zero-crossings *)
+    mutable horizon : float; (* the next horizon *)
+    mutable major : bool; (* integration iff [major = false] *)
+  }
+
+(* A hybrid node is a node that is parameterised by a continuous state *)
+(* all instances points to this global parameter and read/write on it *)
+type ('a, 'b) hnode = cstate -> (time * 'a, 'b) node
+
+type 'b hsimu =
+    Hsim:
+      { alloc : unit -> 's;
+        (* allocate the initial state *)
+        maxsize : 's -> int * int;
+        (* returns the max length of the *)
+        (* cvector and zvector *)
+        csize : 's -> int;
+        (* returns the current length of the continuous state vector *)
+        zsize : 's -> int;
+        (* returns the current length of the zero-crossing vector *)
+        step : 's -> cvec -> dvec -> zinvec -> time -> 'b;
+        (* computes a step *)
+        derivative : 's -> cvec -> dvec -> zinvec -> zoutvec -> time -> unit;
+        (* computes the derivative *)
+        crossings : 's -> cvec -> zinvec -> zoutvec -> time -> unit;
+        (* computes the zero-crossings *)
+        reset : 's -> unit;
+        (* resets the state *)
+        horizon : 's -> time;
+        (* gives the next time horizon *)
+      } -> 'b hsimu
+
+(* a function with type 'a -C-> 'b, when given to a solver, is *)
+(* represented by an OCaml value of type ('a, 'b) hsnode *)
+type ('a, 'b) hsnode =
+    Hnode:
+      { state : 's;
+        (* the discrete state *)
+        zsize : int;
+        (* the maximum size of the zero-crossing vector *)
+        csize : int;
+        (* the maximum size of the continuous state vector (positions) *)
+        derivative : 's -> 'a -> time -> cvec -> dvec -> unit;
+        (* computes the derivative *)
+        crossing : 's -> 'a -> time -> cvec -> zoutvec -> unit;
+        (* computes the derivative *)
+        output : 's -> 'a -> cvec -> 'b;
+        (* computes the zero-crossings *)
+        setroots : 's -> 'a -> cvec -> zinvec -> unit;
+        (* returns the zero-crossings *)
+        majorstep : 's -> time -> cvec -> 'a -> 'b;
+        (* computes a step *)
+        reset : 's -> unit;
+        (* resets the state *)
+        horizon : 's -> time;
+        (* gives the next time horizon *)
+      } -> ('a, 'b) hsnode
+
+ (* An idea suggested by Adrien Guatto, 26/04/2021 *)
+ (* provide a means to the type for input/outputs of nodes *)
+ (* express them with GADT to ensure type safety *)
+ (* type ('a, 'b) node =
+  |  Fun : { step : 'a -> 'b;
+             typ_arg: 'a typ;
+             typ_return: 'b typ
+           } -> ('a, 'b) node
+  | Node :
+      { state : 's; step : 's -> 'a -> 'b * 's;
+        typ_arg: 'a typ;
+        typ_state : 's typ;
+        typ_return: 'b typ } -> ('a, 'b) node
+
+and 'a typ =
+  | Tunit : unit typ
+  | Tarrow : 'a typ * 'b typ -> ('a * 'b) typ
+  | Tint : int -> int typ
+  | Ttuple : 'a typlist -> 'a typ
+  | Tnode : 'a typ * 'b typ -> ('a,'b) node typ
+
+and 'a typlist =
+  | Tnil : unit typlist
+  | Tpair : 'a typ * 'b typlist -> ('a * 'b) typlist
+
+Q1: do it for records? sum types ? How?
+Q2: provide a "type_of" function for every introduced type?
+*)
+
diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml
index 346b1db..d6dd628 100644
--- a/src/lib/hsim/sim.ml
+++ b/src/lib/hsim/sim.ml
@@ -3,202 +3,195 @@ 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 zset csize zsize jump reset =
+      let ms, ss, zin = get_mstate s, get_sstate s, get_zin s in
+      let ms = match zin with Some z -> zset ms z | None -> ms 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 s = set_zin None s in
+        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
+      Utils.dot o, s
+
+    let step_continuous s step cset fout hor =
+      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 stop = min stop (hor ms) in
+      let (h, f, z), ss = step ss (min stop (hor ms)) in
+      let ms = cset ms (f h) in
+      let fy t = f (now +. t) in
+      let fms t = cset ms (fy t) in
+      let fout t = fout ms (i.u (now +. t)) (fy t) in
+      let s, c = match z with
+      | None ->
+          if h >= stop
+          then set_running ~mode:Discrete ~now:h s, Discontinuous
+          else set_running ~now: h s, Continuous
+      | Some _ -> set_running ~mode:Discrete ~now:h s, Discontinuous in
+      let h = h -. now in
+      { h; u=fout; c }, update ms ss (set_zin z s), { h; c; u=fms }
+
+    (** 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
+      (HNode m : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
+      (DNode s : ('y, 'yder, 'zin, 'zout) solver)
+    : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim
+    = let state = get_init m.state s.state in
+      let step_discrete st =
+        let o, s = step_discrete st m.step m.horizon m.fder m.fzer m.cget m.zset
+          m.csize m.zsize m.jump s.reset in
+        Some o, s in
+      let step_continuous st =
+        let o, s, _ = step_continuous st s.step m.cset m.fout m.horizon in
+        Some o, s in
 
-      let step s i =
-        let ms, ss = get_mstate s, get_sstate s in
-        match i, is_running s with
-        | Some i, _ ->
-            let mode, now, stop = Discrete, 0.0, i.length 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 now t) in
-                    let fzer t = model.fzer ms (Utils.offset i 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      = { start=i.start +. now; length=i.length -. now;
-                                   u=Utils.offset i now } in
-                    let mode, stop, now = Continuous, i.length, 0.0 in
-                    update ms ss (set_running ~mode ~input:i ~stop ~now s)
-                  end else set_running ~mode:Continuous s in
-                Some { start=i.start +. now; length=0.0; u=fun _ -> o }, s
+      let 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 -> step_continuous st
+            else None, st 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 }
+
+    let rec run_assert :
+      'a 'b. ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a ->
+        (unit -> ('y, 'yder, 'zin, 'zout) solver) ->
+        ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim =
+    fun (HNodeA m) get_s ->
+      let DNode s = get_s () in
+      let al = List.map (fun a -> run_assert a get_s) m.assertions in
+      let state = get_init m.body.state s.state, al in
+
+      let step_discrete (st, al) =
+        let m=m.body in
+        let o, st =
+          step_discrete st m.step m.horizon m.fder m.fzer m.cget m.zset m.csize
+          m.zsize m.jump s.reset in
+        let al = List.map (fun (DNode a) ->
+          let _, state = a.step a.state @@ Some (Utils.dot @@ get_mstate st) in
+          DNode { a with state }) al in
+        Some o, (st, al) in
+
+      let step_continuous (st, al) =
+        let ({ h; _ } as o), st, u =
+          step_continuous st s.step m.body.cset m.body.fout m.body.horizon in
+        let al = List.map (fun (DNode a) ->
+          (* Step assertions repeatedly until they reach the horizon. *)
+          let rec step s =
+            let o, s = a.step s None in
+            match o with None -> s | Some _ -> step s in
+          let state = step (snd @@ a.step a.state (Some u)) in
+          DNode { a with state }) al in
+        (* Reset the model's state to the reached horizon. *)
+        let st = set_mstate (u.u h) st in
+        Some o, (st, al) in
+
+      let step (st, al) = function
+        | Some i ->
+            let mode, now, stop = Discrete, 0.0, i.h in
+            step_discrete (set_running ~mode ~input:i ~now ~stop st, al)
+        | None ->
+            if is_running st then match get_mode st with
+            | Discrete   -> step_discrete (st, al)
+            | Continuous -> step_continuous (st, al)
+            else None, (st, al) in
+
+      let reset (pm, ps) (st, al) =
+        let ms = m.body.reset pm (get_mstate st) in
+        let ss = s.reset ps (get_sstate st) in
+        update ms ss (set_idle st), al in
+
+      DNode { state; step; reset }
+
+    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) sim
+    = let state = get_init m.state s.state in
+      let step_discrete st =
+        let o, st = step_discrete st m.step m.horizon m.fder m.fzer m.cget
+          m.zset m.csize m.zsize m.jump s.reset in
+        Some o, st in
+      let step_continuous st =
+        let o, st, _ = step_continuous st s.step m.cset m.fout m.horizon in
+        o, st in
+
+      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 (h, f, z), ss = solver.step ss stop in
-                let ms = model.cset ms (f h) in
-                let start = i.start +. now in
-                let fout t = model.fout ms (i.u (now +. t)) (f (now +. t)) in
-                let out = { start; length=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
+              let o, st = step_continuous st in
+              match o.c with
+              | Discontinuous -> Some o, st
+              | Continuous ->
+                  let o', st = step st None in
+                  Some (Utils.concat [o; Option.get o']), st
+            else None, st in
 
-      let reset (pm, ps) s =
-        let ms = model.reset pm (get_mstate s) in
-        let ss = solver.reset ps (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 state = sim.step sim.state (Some input) in
-      let state = match state with None, s -> s | _ -> assert false in
+    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 = sim.step state None in
+        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 run_on_n (DNode n) inputs use =
       ignore @@ List.fold_left (fun state i ->
-        let state = match sim.step state (Some i) with
-        | None, s -> s | _ -> assert false in
+        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 = sim.step state None in
+          let o, state = n.step state None in
           match o with None -> state | Some o -> use o; loop state in
-        loop state) sim.state inputs
+        loop state) n.state inputs
 
-    (** Run the model autonomously until [length], or until the model stops
+    (** Run the model autonomously until [h], or until the model stops
         answering. *)
-    let run_until model solver length =
-      run_on model solver { start = 0.0; length; 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 [length], split in multiple [steps]. *)
-    let run_until_n model solver length steps =
-      let step = length /. (float_of_int steps) in
-      let inputs = List.init steps (fun s ->
-        let start  = float_of_int s *. step in
-        let stop   = min (float_of_int (s+1) *. step) length in
-        { start; length = stop -. start; u = fun _ -> () }) in
-      run_on_n model solver inputs
-  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.length 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 now t) in
-                let fzer t = model.fzer ms (Utils.offset i 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      = { start=i.start +. now; length=i.length -. now;
-                               u=Utils.offset i now } in
-                let mode, stop, now = Continuous, i.length, 0.0 in
-                let s = set_running ~mode ~input:i ~stop ~now s in
-                step (update ms ss s) i
-              else step (set_running ~mode:Continuous s) i in
-            { start = i.start +. now; length = 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 = { start = i.start +. now; length = 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.compose [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
-
-      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
-
-    (** 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))
-
-    (** Run the model autonomously until [length], or until the model stops
-        answering. *)
-    let run_until model solver length =
-      run_on model solver { start = 0.0; length; u = fun _ -> () }
-
-    (** Run the model autonomously until [length], split in multiple [steps]. *)
-    let run_until_n model solver length steps =
-      let step = length /. (float_of_int steps) in
-      let inputs = List.init steps (fun s ->
-        let start  = float_of_int s *. step in
-        let stop   = min (float_of_int (s+1) *. step) length in
-        { start; length = stop -. start; u = fun _ -> () }) in
-      run_on_n model solver inputs
   end
diff --git a/src/lib/hsim/state.ml b/src/lib/hsim/state.ml
index 676e12a..29fd8b2 100644
--- a/src/lib/hsim/state.ml
+++ b/src/lib/hsim/state.ml
@@ -15,59 +15,65 @@ module type SimState =
         - Idle:    waiting for input;
         - Running: currently integrating; in this case, we have access to the
                    step mode, current input, timestamp and stop time. *)
-    type ('a, 'ms, 'ss) state
+    type ('a, 'ms, 'ss, 'zin) state
 
     (** Get the model state. *)
-    val get_mstate : ('a, 'ms, 'ss) state -> 'ms
+    val get_mstate : ('a, 'ms, 'ss, 'zin) state -> 'ms
 
     (** Get the solver state. *)
-    val get_sstate : ('a, 'ms, 'ss) state -> 'ss
+    val get_sstate : ('a, 'ms, 'ss, 'zin) state -> 'ss
+
+    (** Get the last zero-crossing value. *)
+    val get_zin : ('a, 'ms, 'ss, 'zin) state -> 'zin option
 
     (** Get the current step mode.
         ⚠ Should only be called when running (see [is_running]). *)
-    val get_mode : ('a, 'ms, 'ss) state -> mode
+    val get_mode : ('a, 'ms, 'ss, 'zin) state -> mode
 
     (** Get the current input.
         ⚠ Should only be called when running (see [is_running]). *)
-    val get_input : ('a, 'ms, 'ss) state -> 'a value
+    val get_input : ('a, 'ms, 'ss, 'zin) state -> 'a value
 
     (** Get the current timestamp.
         ⚠ Should only be called when running (see [is_running]). *)
-    val get_now : ('a, 'ms, 'ss) state -> time
+    val get_now : ('a, 'ms, 'ss, 'zin) state -> time
 
     (** Get the current stop time.
         ⚠ Should only be called when running (see [is_running]). *)
-    val get_stop : ('a, 'ms, 'ss) state -> time
+    val get_stop : ('a, 'ms, 'ss, 'zin) state -> time
 
     (** Build an initial state. *)
-    val get_init : 'ms -> 'ss -> ('a, 'ms, 'ss) state
+    val get_init : 'ms -> 'ss -> ('a, 'ms, 'ss, 'zin) state
 
     (** Is the simulation running or idle ? *)
-    val is_running : ('a, 'ms, 'ss) state -> bool
+    val is_running : ('a, 'ms, 'ss, 'zin) state -> bool
 
     (** Update the model state. *)
-    val set_mstate : 'ms -> ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state
+    val set_mstate : 'ms -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
 
     (** Update the solver state. *)
-    val set_sstate : 'ss -> ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state
+    val set_sstate : 'ss -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
+
+    (** Update the zero-crossing value. *)
+    val set_zin : 'zin option -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
 
     (** Update both the solver and model states. *)
-    val update : 'ms -> 'ss -> ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state
+    val update : 'ms -> 'ss -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
 
     (** Update the status to running. *)
     val set_running :
       ?mode:mode -> ?input:'a value -> ?now:time -> ?stop:time ->
-        ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state
+        ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
 
     (** Update the status to idle. *)
-    val set_idle : ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state
+    val set_idle : ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
   end
 
 module type SimStateCopy =
   sig
     include SimState
 
-    val copy : ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state
+    val copy : ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
   end
 
 module FunctionalSimState : SimState =
@@ -88,15 +94,17 @@ module FunctionalSimState : SimState =
 
     (** Internal state of the simulation node: model state, solver state and
         current simulation status. *)
-    type ('a, 'ms, 'ss) state =
+    type ('a, 'ms, 'ss, 'zin) state =
       { status : 'a status;                  (** Current simulation status.   *)
         mstate : 'ms;                        (** Model state.                 *)
-        sstate : 'ss }                       (** Solver state.                *)
+        sstate : 'ss;                        (** Solver state.                *)
+        zin    : 'zin option; }              (** Last zero-crossing vector    *)
 
     exception Not_running
 
     let get_mstate state = state.mstate
     let get_sstate state = state.sstate
+    let get_zin state = state.zin
 
     let is_running state =
       match state.status with Running _ -> true | Idle -> false
@@ -120,6 +128,7 @@ module FunctionalSimState : SimState =
 
     let set_mstate mstate state = { state with mstate }
     let set_sstate sstate state = { state with sstate }
+    let set_zin zin state = { state with zin }
 
     let update mstate sstate state = { state with mstate; sstate }
 
@@ -132,7 +141,7 @@ module FunctionalSimState : SimState =
     let get_stop s =
       match s.status with Running r -> r.stop | Idle -> raise Not_running
 
-    let get_init mstate sstate = { status = Idle; mstate; sstate }
+    let get_init mstate sstate = { status = Idle; mstate; sstate; zin = None }
   end
 
 module InPlaceSimState : SimState =
@@ -146,15 +155,17 @@ module InPlaceSimState : SimState =
             mutable stop  : time;
           } -> 'a status
 
-    type ('a, 'ms, 'ss) state =
+    type ('a, 'ms, 'ss, 'zin) state =
       { mutable status : 'a status;
         mutable mstate : 'ms;
-        mutable sstate : 'ss }
+        mutable sstate : 'ss;
+        mutable zin    : 'zin option }
 
     exception Not_running
 
     let get_mstate state = state.mstate
     let get_sstate state = state.sstate
+    let get_zin state = state.zin
 
     let is_running state =
       match state.status with Running _ -> true | Idle -> false
@@ -179,6 +190,7 @@ module InPlaceSimState : SimState =
 
     let set_mstate mstate state = state.mstate <- mstate; state
     let set_sstate sstate state = state.sstate <- sstate; state
+    let set_zin zin state = state.zin <- zin; state
 
     let update mstate sstate state =
       state.mstate <- mstate; state.sstate <- sstate; state
@@ -192,6 +204,6 @@ module InPlaceSimState : SimState =
     let get_stop s =
       match s.status with Running r -> r.stop | Idle -> raise Not_running
 
-    let get_init mstate sstate = { status = Idle; mstate; sstate }
+    let get_init mstate sstate = { status = Idle; mstate; sstate; zin=None }
   end
 
diff --git a/src/lib/hsim/types.ml b/src/lib/hsim/types.ml
index e452591..23e7b52 100644
--- a/src/lib/hsim/types.ml
+++ b/src/lib/hsim/types.ml
@@ -1,69 +1,77 @@
 
 type time = float
+type continuity = Continuous | Discontinuous
 
 (** Input and output values are functions defined on intervals. *)
 type 'a value =
-  { start : time;
-    length : time;                     (* Relative: [end = start + length].   *)
-    u : time -> 'a }                   (* Defined on [[start, end]].          *)
+  { h : time;
+    u : time -> 'a;                    (* Defined on [[0, h]].                *)
+    c : continuity }                   (* Continuity w.r.t. the next value.   *)
 
 (** A time signal is a sequence of possibly absent α-values
-    [{ start; length; u }] where:
-  - [start] and [length] are positive (possibly null) floating-point numbers;
-  - [u: [0, length] -> α] *)
+    [{ h; u }] where:
+  - [h : R⁺]
+  - [u: [0, h] -> α] *)
 type 'a signal = 'a value option
 
+type ('s, 'p, 'a, 'b) drec =
+  { state : 's;
+    step  : 's -> 'a -> 'b * 's;
+    reset : 'p -> 's -> 's }
+
 (** A discrete node. *)
 type ('p, 'a, 'b) dnode =
-  DNode :
-    { state : 's;
-      step  : 's -> 'a -> 'b * 's;
-      reset : 'p -> 's -> 's;
-    } -> ('p, 'a, 'b) dnode
+  DNode : ('s, 'p, 'a, 'b) drec -> ('p, 'a, 'b) dnode
+
+type ('s, 'p, 'a, 'b) drec_c =
+  { state : 's;
+    step  : 's -> 'a -> 'b * 's;
+    reset : 'p -> 's -> 's;
+    copy  : 's -> 's }
 
 (** A discrete node which supports a state copy. *)
 type ('p, 'a, 'b) dnode_c =
-  DNodeC :
-    { state : 's;
-      step  : 's -> 'a -> 'b * 's;
-      reset : 'p -> 's -> 's;
-      copy  : 's -> 's;
-    } -> ('p, 'a, 'b) dnode_c
+  DNodeC : ('s, 'p, 'a, 'b) drec_c -> ('p, 'a, 'b) dnode_c
 
-(** A continuous node. *)
-type ('a, 'b, 'y, 'yder) cnode =
-  CNode :
-    { lsty : 'y;
-      fder : 'a -> 'y -> 'yder;
-      fout : 'a -> 'y -> 'b;
-    } -> ('a, 'b, 'y, 'yder) cnode
+type ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec =
+  { state: 's;
+    step : 's -> 'a -> 'b * 's;        (** Discrete step function.            *)
+    fder : 's -> 'a -> 'y -> 'yder;    (** Continuous derivative function.    *)
+    fout : 's -> 'a -> 'y -> 'b;       (** Continuous output function.        *)
+    fzer : 's -> 'a -> 'y -> 'zout;    (** Continuous zero-crossing function. *)
+    reset : 'p -> 's -> 's;            (** Reset function.                    *)
+    horizon : 's -> time;              (** Next integration horizon.          *)
+    jump : 's -> bool;                 (** Discontinuity flag.                *)
+    cget : 's -> 'y;                   (** Get continuous state.              *)
+    cset : 's -> 'y -> 's;             (** Set continuous state.              *)
+    zset : 's -> 'zin -> 's;           (** Set zero-crossing state.           *)
+    csize : int;
+    zsize : int }
 
 (** A hybrid node. *)
 type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
-  HNode :
-    { state: 's;
-      step : 's -> 'a -> 'b * 's;    (** Discrete step function.            *)
-      fder : 's -> 'a -> 'y -> 'yder; (** Continuous derivative function.    *)
-      fout : 's -> 'a -> 'y -> 'b;    (** Continuous output function.        *)
-      fzer : 's -> 'a -> 'y -> 'zout; (** Continuous zero-crossing function. *)
-      reset : 'p -> 's -> 's;        (** Reset function.                    *)
-      horizon : 's -> time;           (** Next integration horizon.          *)
-      jump : 's -> bool;              (** Discontinuity flag.                *)
-      cget : 's -> 'y;                (** Get continuous state.              *)
-      cset : 's -> 'y -> 's;         (** Set continuous state.              *)
-      zset : 's -> 'zin -> 's;       (** Set zero-crossing state.           *)
-      csize : int;
-      zsize : int;
-    } -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
+  HNode : ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec ->
+          ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
+
+(** A hybrid node with assertions. *)
+type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a =
+  HNodeA : {
+    body : ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec;
+    assertions : ('p, 's, unit, 'y, 'yder, 'zin, 'zout) hnode_a list
+  } -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a
 
 (** The simulation of a hybrid system is a synchronous function on streams of
     functions. *)
-type ('p, 'a, 'b) lazy_sim = ('p, 'a signal, 'b signal) dnode
+type ('p, 'a, 'b) 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
-
-(** Utils *)
+(* Utils *)
 
+(** Consider a node with state copying as a node without state copying. *)
 let d_of_dc (DNodeC { state; step; reset; _ }) = DNode { state; step; reset }
+
+(** Consider a model without assertions as a model with an empty list of
+    assertions. *)
+let a_of_h (HNode body) = HNodeA { body; assertions=[] }
+
+(** Consider a model without its assertions. *)
+let h_of_a (HNodeA { body; _ }) = HNode body
diff --git a/src/lib/hsim/utils.ml b/src/lib/hsim/utils.ml
index e4da67f..76313a1 100644
--- a/src/lib/hsim/utils.ml
+++ b/src/lib/hsim/utils.ml
@@ -1,24 +1,80 @@
 
 open Types
 
-(** Offset the [input.u] function by [now]. *)
-let offset (input : 'a value) (now : time) : time -> 'a =
-  fun t -> input.u ((now -. input.start) +. t)
+let dot v = { h=0.0; c=Discontinuous; u=fun _ -> v }
 
-(**
-Concatenate functions. [
-^                    ^
-|     ---,           |     ---,
-|    ___  `---    =  |    _    `---
-| --'                | --'
-+-------------->     +-------------->]
-*)
-let rec compose = function
+(** Offset the [u] function by [now]. *)
+let offset (u : time -> 'a) (now : time) : time -> 'a =
+  fun t -> u (t +. now)
+
+(** Concatenate functions. *)
+let rec concat = function
   | [] -> raise (Invalid_argument "Cannot concatenate an empty value list")
   | [f] -> f
-  | { start; u; _ } :: l ->
-      let { start=sr; length=lr; u=ur } = compose l in
-      let sw = sr -. start in
-      let length = sw +. lr  in
-      { start; length; u=fun t -> if t < sw then u t else ur (t -. sw) }
+  | { 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) }
 
+(** Sample a function at [n] equidistant points. *)
+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 }
+
+(** Compose two simulations. *)
+let compose_sim
+  (DNode m : ('p, 'a, 'b) sim)
+  (DNode n : ('q, 'b, 'c) 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 }
diff --git a/src/lib/solvers/illinois.ml b/src/lib/solvers/illinois.ml
index 434be45..7e280ce 100644
--- a/src/lib/solvers/illinois.ml
+++ b/src/lib/solvers/illinois.ml
@@ -3,7 +3,9 @@
 (* part of the Zelus standard library.                                        *)
 (* It is implemented with in-place modification of arrays.                    *)
 
-let debug = ref false
+let debug () =
+  (* false *)
+  !Common.Debug.debug
 		
 let printf x = Format.printf x
 
@@ -121,7 +123,7 @@ type t = {
 let reinitialize ({ g; f1 = f1; t1 = t1; _ } as s) t c =
   s.t1 <- t;
   g t1 c f1;   (* fill f1, because it is immediately copied into f0 by next_mesh *)
-  if !debug then (printf "z|---------- init(%.24e, ... ----------@." t;
+  if debug () then (printf "z|---------- init(%.24e, ... ----------@." t;
                   log_limit s.f1);
   s.bothf_valid  <- false
 
@@ -162,7 +164,7 @@ let step ({ g; f0 = f0; f1 = f1; t1 = t1; _ } as s) t c =
   g t c s.f1;
   s.bothf_valid <- true;
 
-  if !debug then
+  if debug () then
     (printf "z|---------- step(%.24e, %.24e)----------@." s.t0 s.t1;
      log_limits s.f0 s.f1)
 
@@ -212,7 +214,7 @@ let find ({ g = g; bothf_valid = bothf_valid;
     dky t_right 0;      (* c = dky_0(t_right); update state *)
     ignore (update_roots calc_zc f_left (get_f_right f_right') roots);
 
-    if !debug then
+    if debug () then
       (printf
           "z|---------- stall(%.24e, %.24e) {interval < %.24e !}--@."
           t_left t_right ttol;
@@ -280,20 +282,20 @@ let find ({ g = g; bothf_valid = bothf_valid;
 
         match check_interval calc_zc f_left f_mid with
         | SearchLeft  ->
-            if !debug then printf "z| (%.24e -- %.24e]   %.24e@."
+            if debug () then printf "z| (%.24e -- %.24e]   %.24e@."
                            t_left t_mid t_right;
             let alpha = if i >= 1 then alpha *. 0.5 else alpha in
             let n_mid = f_mid_from_f_right f_right' in
             seek (t_left, f_left, n_mid,  t_mid,   Some f_mid, alpha, i + 1)
 
         | SearchRight ->
-            if !debug then printf "z|  %.24e   (%.24e -- %.24e]@."
+            if debug () then printf "z|  %.24e   (%.24e -- %.24e]@."
                            t_left t_mid t_right;
             let alpha = if i >= 1 then alpha *. 2.0 else alpha in
             seek (t_mid,  f_mid,  f_left, t_right, f_right',   alpha, i + 1)
 
         | FoundMid    ->
-            if !debug then printf "z|  %.24e   [%.24e]   %.24e@."
+            if debug () then printf "z|  %.24e   [%.24e]   %.24e@."
                            t_left t_mid t_right;
             ignore (update_roots calc_zc f_left f_mid roots);
             let f_tmp = f_mid_from_f_right f_right' in
@@ -303,7 +305,7 @@ let find ({ g = g; bothf_valid = bothf_valid;
 
   if not bothf_valid then (clear_roots roots; assert false)
   else begin
-    if !debug then
+    if debug () then
       printf "z|\nz|---------- find(%.24e, %.24e)----------@." t0 t1;
 
     match check_interval calc_zc f0 f1 with
@@ -314,7 +316,7 @@ let find ({ g = g; bothf_valid = bothf_valid;
       end
 
     | FoundMid    -> begin
-        if !debug then printf "z| zero-crossing at limit (%.24e)@." t1;
+        if debug () then printf "z| zero-crossing at limit (%.24e)@." t1;
         ignore (update_roots calc_zc f0 f1 roots);
         s.bothf_valid <- false;
         t1
diff --git a/src/lib/solvers/odexx.ml b/src/lib/solvers/odexx.ml
index 31637aa..30252be 100644
--- a/src/lib/solvers/odexx.ml
+++ b/src/lib/solvers/odexx.ml
@@ -51,7 +51,9 @@ module GenericODE (Butcher : BUTCHER_TABLEAU) : STATE_ODE_SOLVER =
 struct (* {{{1 *)
   open Bigarray
 
-  let debug = ref false (* !Debug.debug *)
+  let debug () =
+    false
+    (* !Common.Debug.debug *)
 
   let pow = 1.0 /. float(Butcher.order)
 
@@ -274,7 +276,7 @@ struct (* {{{1 *)
            "odexx: step size < min step size (\n       now=%.24e\n         h=%.24e\n< min_step=%.24e)"
            t h s.min_step);
 
-    if !debug then Printf.printf "s|\ns|----------step(%.24e)----------\n" max_t;
+   if debug () then Printf.printf "s|\ns|----------step(%.24e)----------\n" max_t;
 
     let rec onestep (alreadyfailed: bool) h =
 
@@ -288,11 +290,11 @@ struct (* {{{1 *)
       let tnew = if finished then max_t else t +. h *. (mA maxK) in
       mapinto ynew (make_newval y k maxK);
       f tnew ynew k.(maxK);
-      if !debug then log_step t y k.(0) tnew ynew k.(maxK);
+     if debug () then log_step t y k.(0) tnew ynew k.(maxK);
 
       let err = h *. calculate_error (abs_tol /. rel_tol) k y ynew in
       if err > rel_tol then begin
-        if !debug then Printf.printf "s| error exceeds tolerance\n";
+       if debug () then Printf.printf "s| error exceeds tolerance\n";
 
         if h <= hmin then failwith
             (Printf.sprintf "Error (%e) > relative tolerance (%e) at t=%e"
diff --git a/src/lib/solvers/statefulZ.ml b/src/lib/solvers/statefulZ.ml
index 2b65544..b642a08 100644
--- a/src/lib/solvers/statefulZ.ml
+++ b/src/lib/solvers/statefulZ.ml
@@ -15,7 +15,8 @@ module Functional =
           vec = zmake 0 } in
 
       let reset { fzer; init; size } { vec; _ } =
-        let fzer t cvec zout = let zout' = fzer t cvec in blit zout' zout in
+        let fzer t cvec zout =
+          let zout' = fzer t cvec in blit zout' zout in
         { state = initialize size fzer init;
           vec = if length vec = size then vec else zmake size } in
 
@@ -58,7 +59,9 @@ module InPlace =
           (h, Some vec), s
         else (h, None), s in
 
-      let copy _ = raise Common.Errors.TODO in
+      let copy s =
+        let vec = zmake (length s.vec) in
+        blit s.vec vec; s.vec <- vec; s in
 
       DNodeC { state; step; reset; copy }
   end