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 fa0b780..ecc9b8f 100644
--- a/doc/notes.typ
+++ b/doc/notes.typ
@@ -335,7 +335,7 @@ composition breaks this:
   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                     & |
@@ -353,7 +353,7 @@ needed. This is done through the following composition:
   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 & |
@@ -620,11 +620,11 @@ simulation loop as follows:
   models has much in common with code generation for synchronous languages.
 ])
 
-In [Branicky, 1995b], dynamical systems are defined as follows: 
+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 
+  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$,
@@ -694,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..72e804e 100644
--- a/exm/ball.ml
+++ b/exm/ball.ml
@@ -40,7 +40,7 @@ let step ({ zin; lx; _ } as s) zfalse =
     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 }
 
-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..4eecb1e 100644
--- a/exm/dune
+++ b/exm/dune
@@ -1,3 +1,5 @@
 (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.ml b/exm/zelus/ballz/ballz.ml
new file mode 100644
index 0000000..08df6f6
--- /dev/null
+++ b/exm/zelus/ballz/ballz.ml
@@ -0,0 +1,131 @@
+(* The Zelus compiler, version 2.2-dev
+  (2025-06-16-15:24) *)
+open Common
+open Ztypes
+open Solvers
+
+let (+=) r v = r := !r + v
+
+let g = 9.81
+
+let y0 = 50.
+
+let y'0 = 0.
+
+type ball =
+  { mutable major : bool;
+    mutable h : float;
+    mutable init : bool;
+    mutable z : (bool, float) zerocrossing;
+    mutable y' : float continuous;
+    mutable y : float continuous }
+
+let ball cstate =
+
+  let ball_alloc _ =
+    cstate.cmax <- cstate.cmax + 2;
+    cstate.zmax <- cstate.zmax + 1;
+    { major = false;
+      h = 42.;
+      init = false;
+      z = { zin = false; zout = 1. };
+      y' = { pos = 42.; der = 0. };
+      y  = { pos = 42.; der = 0. }
+    } in
+
+  let ball_step self (_, ()) =
+    let cidx = cstate.cindex in let cpos = ref cidx in
+    let zidx = cstate.zindex in let zpos = ref zidx in
+    cstate.cindex <- cstate.cindex + 2;
+    cstate.zindex <- cstate.zindex + 1;
+    self.major <- cstate.major;
+    if cstate.major then
+      for i = cidx to 1 do Zls.set cstate.dvec i 0. done
+    else begin
+      self.y'.pos <- Zls.get cstate.cvec !cpos; cpos += 1;
+      self.y.pos <- Zls.get cstate.cvec !cpos; cpos += 1
+    end;
+    let res0, res1 =
+      let encore = ref false in
+      if self.init then self.y'.pos <- y'0;
+      let last_y' = self.y'.pos in
+      if self.z.zin then begin
+        encore := true; self.y'.pos <- -0.8 *. last_y'
+      end;
+      self.h <- if !encore then 0. else infinity;
+      if self.init then self.y.pos <- y0;
+      self.init <- false;
+      self.z.zout <- -. self.y.pos;
+      self.y'.der <- -. g;
+      self.y.der <- self.y'.pos;
+      self.y.pos, self.y.der
+    in
+    cstate.horizon <- min cstate.horizon self.h;
+    cpos := cidx;
+    if cstate.major then begin
+      Zls.set cstate.cvec !cpos self.y'.pos; cpos += 1;
+      Zls.set cstate.cvec !cpos self.y.pos; cpos += 1;
+      self.z.zin <- false
+    end else begin
+      self.z.zin <- Zls.get_zin cstate.zinvec !zpos; zpos += 1;
+      zpos := zidx;
+      Zls.set cstate.zoutvec !zpos self.z.zout; zpos += 1;
+      Zls.set cstate.dvec    !cpos self.y'.der; cpos += 1;
+      Zls.set cstate.dvec    !cpos self.y.der;  cpos += 1
+    end;
+    Bigarray.(Array1.of_array Float64 c_layout [| res0; res1 |]) in
+
+  let ball_reset self = self.init <- true in
+
+  Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset }
+
+type ('f , 'e , 'd , 'c , 'b , 'a) _main =
+  { mutable i_73 : 'f ;
+    mutable major_62 : 'e ;
+    mutable h_72 : 'd ;
+    mutable i_70 : 'c ; mutable h_68 : 'b ; mutable t_63 : 'a }
+
+let main (cstate_80:Ztypes.cstate) =
+  let Node { alloc = i_73_alloc; step = i_73_step ; reset = i_73_reset } = ball
+  cstate_80 in
+  let main_alloc _ =
+    cstate_80.cmax <- (+) cstate_80.cmax  1;
+    { major_62 = false ;
+      h_72 = 42. ;
+      i_70 = (false:bool) ;
+      h_68 = (42.:float) ; t_63 = { pos = 42.; der = 0. };
+      i_73 = i_73_alloc () (* continuous *)  } in
+  let main_step self ((time_61:float) , ()) =
+    ((let (cindex_81:int) = cstate_80.cindex in
+      let cpos_83 = ref (cindex_81:int) in
+      cstate_80.cindex <- (+) cstate_80.cindex  1 ;
+      self.major_62 <- cstate_80.major ;
+      (if cstate_80.major then
+       for i_1 = cindex_81 to 0 do Zls.set cstate_80.dvec  i_1  0. done
+       else ((self.t_63.pos <- Zls.get cstate_80.cvec  !cpos_83 ;
+              cpos_83 := (+) !cpos_83  1))) ;
+      (let (result_85) =
+           let h_71 = ref (infinity:float) in
+           (if self.i_70 then self.h_68 <- (+.) time_61  0.) ;
+           (let (z_69:bool) = (&&) self.major_62  ((>=) time_61  self.h_68) in
+            self.h_68 <- (if z_69 then (+.) self.h_68  0.01 else self.h_68) ;
+            h_71 := min !h_71  self.h_68 ;
+            self.h_72 <- !h_71 ;
+            self.i_70 <- false ;
+            self.t_63.der <- 1. ;
+            (let (y_64:float) = (i_73_step self.i_73 (time_61 , ())).{0} in
+             (begin match z_69 with
+                    | true -> Printf.printf "%.10e\t%.10e\n" self.t_63.pos y_64
+                    | _ -> () end) ;
+             Bigarray.(Array1.create Float64 c_layout 0))) in
+       cstate_80.horizon <- min cstate_80.horizon  self.h_72 ;
+       cpos_83 := cindex_81 ;
+       (if cstate_80.major then
+        (((Zls.set cstate_80.cvec  !cpos_83  self.t_63.pos ;
+           cpos_83 := (+) !cpos_83  1)))
+        else (((Zls.set cstate_80.dvec  !cpos_83  self.t_63.der ;
+                cpos_83 := (+) !cpos_83  1)))) ; result_85))) in
+  let main_reset self  =
+    ((self.i_70 <- true ; self.t_63.pos <- 0. ; i_73_reset self.i_73 ):
+    unit) in
+  Node { alloc = main_alloc; step = main_step ; reset = main_reset }
diff --git a/exm/zelus/ballz/ballz.ml.bak b/exm/zelus/ballz/ballz.ml.bak
new file mode 100644
index 0000000..d248fc0
--- /dev/null
+++ b/exm/zelus/ballz/ballz.ml.bak
@@ -0,0 +1,137 @@
+(* The Zelus compiler, version 2.2-dev
+  (2025-06-16-15:24) *)
+open Common
+open Ztypes
+open Solvers
+
+let g = 9.81
+
+let y0 = 50.
+
+let y'0 = 0.
+
+type ('g , 'f , 'e , 'd , 'c , 'b , 'a) _ball =
+  { mutable major_50 : 'g ;
+    mutable h_60 : 'f ;
+    mutable h_58 : 'e ;
+    mutable i_56 : 'd ;
+    mutable x_55 : 'c ; mutable y'_52 : 'b ; mutable y_51 : 'a }
+
+let ball (cstate_74:Ztypes.cstate) =
+
+  let ball_alloc _ =
+    cstate_74.cmax <- (+) cstate_74.cmax  2 ;
+    cstate_74.zmax <- (+) cstate_74.zmax  1;
+    { major_50 = false ;
+      h_60 = 42. ;
+      h_58 = (42.:float) ;
+      i_56 = (false:bool) ;
+      x_55 = { zin = false; zout = 1. } ;
+      y'_52 = { pos = 42.; der = 0. } ; y_51 = { pos = 42.; der = 0. } } in
+  let ball_step self ((_time_49:float) , ((y0_48:float) , (y'0_47:float))) =
+    Printf.printf "STEP (%d)\n" cstate_74.cindex;
+    ((let (cindex_75:int) = cstate_74.cindex in
+      let cpos_77 = ref (cindex_75:int) in
+      let (zindex_76:int) = cstate_74.zindex in
+      let zpos_78 = ref (zindex_76:int) in
+      cstate_74.cindex <- (+) cstate_74.cindex  2 ;
+      cstate_74.zindex <- (+) cstate_74.zindex  1 ;
+      self.major_50 <- cstate_74.major ;
+      (if cstate_74.major then
+       for i_1 = cindex_75 to 1 do Zls.set cstate_74.dvec  i_1  0. done
+       else ((self.y'_52.pos <- Zls.get cstate_74.cvec  !cpos_77 ;
+              cpos_77 := (+) !cpos_77  1) ;
+             (self.y_51.pos <- Zls.get cstate_74.cvec  !cpos_77 ;
+              cpos_77 := (+) !cpos_77  1))) ;
+      (let (result_79:float) =
+           let h_59 = ref (infinity:float) in
+           let encore_57 = ref (false:bool) in
+           (if self.i_56 then self.y'_52.pos <- y'0_47) ;
+           (let (l_54:float) = self.y'_52.pos in
+            (begin match self.x_55.zin with
+                   | true ->
+                       encore_57 := true ;
+                       self.y'_52.pos <- ( *. ) (-0.8)  l_54 | _ -> ()  end)
+            ;
+            self.h_58 <- (if !encore_57 then 0. else infinity) ;
+            h_59 := min !h_59  self.h_58 ;
+            self.h_60 <- !h_59 ;
+            (if self.i_56 then self.y_51.pos <- y0_48) ;
+            self.i_56 <- false ;
+            self.x_55.zout <- (~-.) self.y_51.pos ;
+            self.y'_52.der <- (~-.) g ;
+            self.y_51.der <- self.y'_52.pos ; self.y_51.pos) in
+       cstate_74.horizon <- min cstate_74.horizon  self.h_60 ;
+       cpos_77 := cindex_75 ;
+       (if cstate_74.major then
+        (((Printf.printf "idx: %d\n" !cpos_77;
+           Zls.set cstate_74.cvec  !cpos_77  self.y'_52.pos ;
+           cpos_77 := (+) !cpos_77  1) ;
+          (Zls.set cstate_74.cvec  !cpos_77  self.y_51.pos ;
+           cpos_77 := (+) !cpos_77  1)) ; ((self.x_55.zin <- false)))
+        else (((self.x_55.zin <- Zls.get_zin cstate_74.zinvec  !zpos_78 ;
+                zpos_78 := (+) !zpos_78  1)) ;
+              zpos_78 := zindex_76 ;
+              ((Zls.set cstate_74.zoutvec  !zpos_78  self.x_55.zout ;
+                zpos_78 := (+) !zpos_78  1)) ;
+              ((Zls.set cstate_74.dvec  !cpos_77  self.y'_52.der ;
+                cpos_77 := (+) !cpos_77  1) ;
+               (Zls.set cstate_74.dvec  !cpos_77  self.y_51.der ;
+                cpos_77 := (+) !cpos_77  1)))) ; result_79)):float) in
+  let ball_reset self  =
+    (self.i_56 <- true:unit) in
+  Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset }
+type ('f , 'e , 'd , 'c , 'b , 'a) _main =
+  { mutable main_ball : 'f ;
+    mutable main_major : 'e ;
+    mutable h_72 : 'd; mutable i_70 : 'c; mutable h : 'b;
+    mutable t_63 : 'a }
+
+let main cs =
+  let Node
+    { alloc = ball_alloc;
+      step = ball_step;
+      reset = ball_reset } = ball cs in
+
+  let main_alloc _ =
+    cs.cmax <- cs.cmax + 1;
+    { main_major = false;
+      h_72 = 42.0; i_70 = false; h = 42.0;
+      t_63 = { pos = 42.; der = 0. };
+      main_ball = ball_alloc () } in
+
+  let main_step self (time, ()) =
+    let cindex = cs.cindex in
+    let cpos = ref cindex in
+    Printf.printf "main:cindex: %d\n" cs.cindex;
+    cs.cindex <- cs.cindex + 1;
+    self.main_major <- cs.major;
+    if cs.major then for i = cindex to 0 do Zls.set cs.dvec i 0. done
+    else begin self.t_63.pos <- Zls.get cs.cvec !cpos; cpos := !cpos + 1 end;
+    let result =
+      if self.i_70 then self.h <- time;
+      let z = self.main_major && (time >= self.h) in
+      if z then self.h <- self.h +. 0.01;
+      self.h_72 <- min infinity self.h;
+      self.i_70 <- false;
+      self.t_63.der <- 1.;
+      let y_64 = ball_step self.main_ball (time, (y0, y'0)) in
+      if z then begin
+        print_float self.t_63.pos;
+        print_string "\t";
+        print_float y_64;
+        print_newline ()
+      end;
+      Bigarray.(Array1.create Float64 c_layout 0) in
+    cs.horizon <- min cs.horizon self.h_72;
+    cpos := cindex;
+    if cs.major then Zls.set cs.cvec !cpos self.t_63.pos
+                else Zls.set cs.dvec !cpos self.t_63.der;
+    result in
+
+  let main_reset self =
+    self.i_70 <- true;
+    self.t_63.pos <- 0.;
+    ball_reset self.main_ball in
+
+  Node { alloc = main_alloc; step = main_step ; reset = main_reset }
diff --git a/exm/zelus/ballz/ballz.ml.bak2 b/exm/zelus/ballz/ballz.ml.bak2
new file mode 100644
index 0000000..05842bb
--- /dev/null
+++ b/exm/zelus/ballz/ballz.ml.bak2
@@ -0,0 +1,135 @@
+(* The Zelus compiler, version 2.2-dev
+  (2025-06-16-15:24) *)
+open Common
+open Ztypes
+open Solvers
+
+let g = 9.81
+
+let y0 = 50.
+
+let y'0 = 0.
+
+type ('g , 'f , 'e , 'd , 'c , 'b , 'a) _ball =
+  { mutable major_50 : 'g ;
+    mutable h_60 : 'f ;
+    mutable h_58 : 'e ;
+    mutable i_56 : 'd ;
+    mutable x_55 : 'c ; mutable y'_52 : 'b ; mutable y_51 : 'a }
+
+let ball (cstate_74:Ztypes.cstate) =
+
+  let ball_alloc _ =
+    cstate_74.cmax <- (+) cstate_74.cmax  2 ;
+    cstate_74.zmax <- (+) cstate_74.zmax  1;
+    { major_50 = false ;
+      h_60 = 42. ;
+      h_58 = (42.:float) ;
+      i_56 = (false:bool) ;
+      x_55 = { zin = false; zout = 1. } ;
+      y'_52 = { pos = 42.; der = 0. } ; y_51 = { pos = 42.; der = 0. } } in
+
+  let ball_step self ((_time_49:float) , ()) =
+    ((let (cindex_75:int) = cstate_74.cindex in
+      let cpos_77 = ref (cindex_75:int) in
+      let (zindex_76:int) = cstate_74.zindex in
+      let zpos_78 = ref (zindex_76:int) in
+      cstate_74.cindex <- (+) cstate_74.cindex  2 ;
+      cstate_74.zindex <- (+) cstate_74.zindex  1 ;
+      self.major_50 <- cstate_74.major ;
+      (if cstate_74.major then
+       for i_1 = cindex_75 to 1 do Zls.set cstate_74.dvec  i_1  0. done
+       else ((self.y'_52.pos <- Zls.get cstate_74.cvec  !cpos_77 ;
+              cpos_77 := (+) !cpos_77  1) ;
+             (self.y_51.pos <- Zls.get cstate_74.cvec  !cpos_77 ;
+              cpos_77 := (+) !cpos_77  1))) ;
+      (let (result_79:float) =
+           let h_59 = ref (infinity:float) in
+           let encore_57 = ref (false:bool) in
+           (if self.i_56 then self.y'_52.pos <- y'0) ;
+           (let (l_54:float) = self.y'_52.pos in
+            begin match self.x_55.zin with
+                   | true ->
+                       encore_57 := true ;
+                       self.y'_52.pos <- ( *. ) (-0.8)  l_54 | _ -> ()  end;
+            self.h_58 <- (if !encore_57 then 0. else infinity) ;
+            h_59 := min !h_59  self.h_58 ;
+            self.h_60 <- !h_59 ;
+            (if self.i_56 then self.y_51.pos <- y0) ;
+            self.i_56 <- false ;
+            self.x_55.zout <- (~-.) self.y_51.pos ;
+            self.y'_52.der <- (~-.) g ;
+            self.y_51.der <- self.y'_52.pos ; self.y_51.pos) in
+       cstate_74.horizon <- min cstate_74.horizon  self.h_60 ;
+       cpos_77 := cindex_75 ;
+       (if cstate_74.major then
+        (((Zls.set cstate_74.cvec  !cpos_77  self.y'_52.pos ;
+           cpos_77 := (+) !cpos_77  1) ;
+          (Zls.set cstate_74.cvec  !cpos_77  self.y_51.pos ;
+           cpos_77 := (+) !cpos_77  1)) ; ((self.x_55.zin <- false)))
+        else (((self.x_55.zin <- Zls.get_zin cstate_74.zinvec  !zpos_78 ;
+                zpos_78 := (+) !zpos_78  1)) ;
+              zpos_78 := zindex_76 ;
+              ((Zls.set cstate_74.zoutvec  !zpos_78  self.x_55.zout ;
+                zpos_78 := (+) !zpos_78  1)) ;
+              ((Zls.set cstate_74.dvec  !cpos_77  self.y'_52.der ;
+                cpos_77 := (+) !cpos_77  1) ;
+               (Zls.set cstate_74.dvec  !cpos_77  self.y_51.der ;
+                cpos_77 := (+) !cpos_77  1)))) ;
+        Bigarray.(Array1.of_array Float64 c_layout [| result_79 |])))) in
+
+  let ball_reset self  =
+    (self.i_56 <- true:unit) in
+
+  Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset }
+
+type ('f , 'e , 'd , 'c , 'b , 'a) _main =
+  { mutable i_73 : 'f ;
+    mutable major_62 : 'e ;
+    mutable h_72 : 'd ;
+    mutable i_70 : 'c ; mutable h_68 : 'b ; mutable t_63 : 'a }
+
+let main (cstate_80:Ztypes.cstate) =
+  let Node { alloc = i_73_alloc; step = i_73_step ; reset = i_73_reset } = ball
+  cstate_80 in
+  let main_alloc _ =
+    cstate_80.cmax <- (+) cstate_80.cmax  1;
+    { major_62 = false ;
+      h_72 = 42. ;
+      i_70 = (false:bool) ;
+      h_68 = (42.:float) ; t_63 = { pos = 42.; der = 0. };
+      i_73 = i_73_alloc () (* continuous *)  } in
+  let main_step self ((time_61:float) , ()) =
+    ((let (cindex_81:int) = cstate_80.cindex in
+      let cpos_83 = ref (cindex_81:int) in
+      cstate_80.cindex <- (+) cstate_80.cindex  1 ;
+      self.major_62 <- cstate_80.major ;
+      (if cstate_80.major then
+       for i_1 = cindex_81 to 0 do Zls.set cstate_80.dvec  i_1  0. done
+       else ((self.t_63.pos <- Zls.get cstate_80.cvec  !cpos_83 ;
+              cpos_83 := (+) !cpos_83  1))) ;
+      (let (result_85) =
+           let h_71 = ref (infinity:float) in
+           (if self.i_70 then self.h_68 <- (+.) time_61  0.) ;
+           (let (z_69:bool) = (&&) self.major_62  ((>=) time_61  self.h_68) in
+            self.h_68 <- (if z_69 then (+.) self.h_68  0.01 else self.h_68) ;
+            h_71 := min !h_71  self.h_68 ;
+            self.h_72 <- !h_71 ;
+            self.i_70 <- false ;
+            self.t_63.der <- 1. ;
+            (let (y_64:float) = (i_73_step self.i_73 (time_61 , ())).{0} in
+             (begin match z_69 with
+                    | true -> Printf.printf "%.10e\t%.10e\n" self.t_63.pos y_64
+                    | _ -> () end) ;
+             Bigarray.(Array1.create Float64 c_layout 0))) in
+       cstate_80.horizon <- min cstate_80.horizon  self.h_72 ;
+       cpos_83 := cindex_81 ;
+       (if cstate_80.major then
+        (((Zls.set cstate_80.cvec  !cpos_83  self.t_63.pos ;
+           cpos_83 := (+) !cpos_83  1)))
+        else (((Zls.set cstate_80.dvec  !cpos_83  self.t_63.der ;
+                cpos_83 := (+) !cpos_83  1)))) ; result_85))) in
+  let main_reset self  =
+    ((self.i_70 <- true ; self.t_63.pos <- 0. ; i_73_reset self.i_73 ):
+    unit) in
+  Node { alloc = main_alloc; step = main_step ; reset = main_reset }
diff --git a/exm/zelus/ballz/ballz.zci b/exm/zelus/ballz/ballz.zci
new file mode 100644
index 0000000..08fa125
Binary files /dev/null and b/exm/zelus/ballz/ballz.zci differ
diff --git a/exm/zelus/ballz/ballz.zls b/exm/zelus/ballz/ballz.zls
new file mode 100644
index 0000000..fbaef0a
--- /dev/null
+++ b/exm/zelus/ballz/ballz.zls
@@ -0,0 +1,19 @@
+let g = 9.81
+let y0 = 50.0
+let y'0 = 0.0
+
+let hybrid ball (y0, y'0) = y 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 y = ball (y0, y'0) in
+  let z = period(0.01) in
+  present z -> (
+    print_float t;
+    print_string "\t";
+    print_float y;
+    print_newline ()
+  ); ()
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/sincosz.ml b/exm/zelus/sincos/sincosz.ml
new file mode 100644
index 0000000..ebbd09d
--- /dev/null
+++ b/exm/zelus/sincos/sincosz.ml
@@ -0,0 +1,45 @@
+(* The Zelus compiler, version 2.2-dev
+  (2025-06-16-15:24) *)
+open Common
+open Ztypes
+open Solvers
+
+type ('c , 'b , 'a) _f =
+  { mutable major_11 : 'c ; mutable sin_13 : 'b ; mutable cos_12 : 'a }
+
+let f (cstate_14:Ztypes.cstate) =
+
+  let f_alloc _ =
+    cstate_14.cmax <- (+) cstate_14.cmax  2;
+    { major_11 = false ;
+      sin_13 = { pos = 42.; der = 0. } ; cos_12 = { pos = 42.; der = 0. } } in
+  let f_step self ((_time_10:float) , ()) =
+    ((let (cindex_15:int) = cstate_14.cindex in
+      let cpos_17 = ref (cindex_15:int) in
+      cstate_14.cindex <- (+) cstate_14.cindex  2 ;
+      self.major_11 <- cstate_14.major ;
+      (if cstate_14.major then
+       for i_1 = cindex_15 to 1 do Zls.set cstate_14.dvec  i_1  0. done
+       else ((self.sin_13.pos <- Zls.get cstate_14.cvec  !cpos_17 ;
+              cpos_17 := (+) !cpos_17  1) ;
+             (self.cos_12.pos <- Zls.get cstate_14.cvec  !cpos_17 ;
+              cpos_17 := (+) !cpos_17  1))) ;
+      (let (result_19) =
+           self.cos_12.der <- (~-.) self.sin_13.pos ;
+           self.sin_13.der <- self.cos_12.pos ;
+           Bigarray.(Array1.of_array Float64 c_layout
+             [| self.sin_13.pos; self.cos_12.pos |]) in
+       cpos_17 := cindex_15 ;
+       (if cstate_14.major then
+        (((Zls.set cstate_14.cvec  !cpos_17  self.sin_13.pos ;
+           cpos_17 := (+) !cpos_17  1) ;
+          (Zls.set cstate_14.cvec  !cpos_17  self.cos_12.pos ;
+           cpos_17 := (+) !cpos_17  1)))
+        else (((Zls.set cstate_14.dvec  !cpos_17  self.sin_13.der ;
+                cpos_17 := (+) !cpos_17  1) ;
+               (Zls.set cstate_14.dvec  !cpos_17  self.cos_12.der ;
+                cpos_17 := (+) !cpos_17  1)))) ; result_19))) in
+
+  let f_reset self  =
+    ((self.sin_13.pos <- 0. ; self.cos_12.pos <- 1.):unit) in
+  Node { alloc = f_alloc; step = f_step ; reset = f_reset }
diff --git a/exm/zelus/sincos/sincosz.zci b/exm/zelus/sincos/sincosz.zci
new file mode 100644
index 0000000..e531551
Binary files /dev/null and b/exm/zelus/sincos/sincosz.zci differ
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/src/bin/lift.ml b/src/bin/lift.ml
new file mode 100644
index 0000000..6b4e36d
--- /dev/null
+++ b/src/bin/lift.ml
@@ -0,0 +1,119 @@
+
+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 =
+    Common.Debug.print "LIFT :: Calling [fzer]";
+    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 044e3ba..de91d4c 100644
--- a/src/bin/main.ml
+++ b/src/bin/main.ml
@@ -4,17 +4,25 @@ open Solvers
 open Examples
 open Common
 open Types
+open Lift
 
-let sample   = ref 10
-let stop     = ref 30.0
+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 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 =
@@ -30,25 +38,43 @@ let opts = [
   "-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";
 ]
 
 let errmsg = "Usage: " ^ Sys.executable_name ^ " [OPTIONS] MODEL\nOptions are:"
 
 let () = try Arg.parse (Arg.align opts) set_model errmsg with _ -> exit 2
 
+let args = List.rev !modelargs
+let () = ignore lift
+
 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 "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 "ballz"  -> lift Ballz.ball
+  | Some "sincosz" -> lift Sincosz.f
   | 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
     let open StatefulSundials in
@@ -57,7 +83,7 @@ let sim =
     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.print !sample (run m s))
+    run_until_n (output !sample (run m s))
   else
     let open StatefulRK45 in
     let c = if !inplace then InPlace.csolve else Functional.csolve in
@@ -66,7 +92,7 @@ let sim =
     let s = Solver.solver_c c z in
     let open Sim.Sim(val st) in
     let n = if !accel then accelerate m s else run m (d_of_dc s) in
-    run_until_n (Output.print !sample n)
+    run_until_n (output !sample n)
 
 let () = sim !stop !steps ignore
 
diff --git a/src/bin/output.ml b/src/bin/output.ml
index c397aa9..e60ea97 100644
--- a/src/bin/output.ml
+++ b/src/bin/output.ml
@@ -1,7 +1,6 @@
 
 open Hsim.Types
 open Hsim.Utils
-open Common
 
 let print_entry y t =
   let n = Bigarray.Array1.dim y in
@@ -13,6 +12,16 @@ let print_entry y t =
   Format.printf "\n";
   flush stdout
 
+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 =
@@ -20,8 +29,16 @@ let print_sample samples ({ h; u; _ }, now) =
     else if i = samples then print_entry (u h) (now +. h)
     else let t = float_of_int i *. step in
     (print_entry (u t) (now +. t); loop (i+1)) in
-  if h <= 0.0 then begin Debug.print "D: "; print_entry (u 0.0) now end
-  else begin Debug.print "C: "; loop 0 end
+  if h <= 0.0 then print_entry (u 0.0) now else loop 0
+
+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
@@ -30,3 +47,18 @@ let print_limits { 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 87406bc..ac8ed4f 100644
--- a/src/lib/hsim/sim.ml
+++ b/src/lib/hsim/sim.ml
@@ -8,6 +8,7 @@ module Sim (S : SimState) =
     include S
 
     let step_discrete s step hor fder fzer cget csize zsize jump reset =
+      Common.Debug.print "SIMU :: DISCRETE :: start";
       let ms, ss = get_mstate s, get_sstate s in
       let i, now, stop = get_input s, get_now s, get_stop s in
       let o, ms = step ms (i.u now) in
@@ -26,9 +27,11 @@ module Sim (S : SimState) =
           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
+      Common.Debug.print "SIMU :: DISCRETE :: end";
       Utils.dot o, s
 
     let step_continuous s step cset fout zset =
+      Common.Debug.print "SIMU :: CONTINUOUS :: start";
       let ms, ss = get_mstate s, get_sstate s in
       let i, now, stop = get_input s, get_now s, get_stop s in
       let (h, f, z), ss = step ss stop in
@@ -46,6 +49,7 @@ module Sim (S : SimState) =
           let s = set_running ~mode:Discrete ~now:h s in
           update (zset ms z) ss s, Discontinuous in
       let h = h -. now in
+      Common.Debug.print "SIMU :: CONTINUOUS :: end";
       { h; u=fout; c }, s, { h; c; u=fms }
 
     (** Simulation of a model with any solver. *)
@@ -55,7 +59,7 @@ module Sim (S : SimState) =
     : ('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 
+        let o, s = step_discrete st m.step m.horizon m.fder m.fzer m.cget
           m.csize m.zsize m.jump s.reset in
         Some o, s in
       let step_continuous st =
@@ -97,7 +101,7 @@ module Sim (S : SimState) =
           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.zset in
diff --git a/src/lib/hsim/types.ml b/src/lib/hsim/types.ml
index 705f20d..23e7b52 100644
--- a/src/lib/hsim/types.ml
+++ b/src/lib/hsim/types.ml
@@ -69,7 +69,7 @@ type ('p, 'a, 'b) sim = ('p, 'a signal, 'b signal) dnode
 (** 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 
+(** Consider a model without assertions as a model with an empty list of
     assertions. *)
 let a_of_h (HNode body) = HNodeA { body; assertions=[] }
 
diff --git a/src/lib/solvers/illinois.ml b/src/lib/solvers/illinois.ml
index 434be45..3f5cdb6 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
 
@@ -152,6 +154,7 @@ let num_roots { f0; _ } = Zls.length f0
 
 (* f0/t0 take the previous values of f1/t1, f1/t1 are refreshed by g *)
 let step ({ g; f0 = f0; f1 = f1; t1 = t1; _ } as s) t c =
+  Common.Debug.print "ZSOL :: Calling [step]";
   (* swap f0 and f1; f0 takes the previous value of f1 *)
   s.f0 <- f1;
   s.t0 <- t1;
@@ -162,7 +165,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 +215,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 +283,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 +306,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 +317,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/statefulRK45.ml b/src/lib/solvers/statefulRK45.ml
index d85658e..997d0ee 100644
--- a/src/lib/solvers/statefulRK45.ml
+++ b/src/lib/solvers/statefulRK45.ml
@@ -22,6 +22,8 @@ module Functional =
         { state; vec = init } in
 
       let step ({ state ; vec=v } as s) h =
+        Common.Debug.print "SOLVER STEP";
+        Common.Debug.print_entry v;
         let y_nv = vec v in
         let h = step state h y_nv in
         let state = copy state in
diff --git a/src/lib/solvers/statefulZ.ml b/src/lib/solvers/statefulZ.ml
index cc8c255..f6dc60b 100644
--- a/src/lib/solvers/statefulZ.ml
+++ b/src/lib/solvers/statefulZ.ml
@@ -15,7 +15,10 @@ 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
+        Common.Debug.print "ZSolver Reset";
+        Common.Debug.print_entry init;
         { state = initialize size fzer init;
           vec = if length vec = size then vec else zmake size } in