From 6cec3d6c5dd13315aacc1cc0c0d0dd674f026229 Mon Sep 17 00:00:00 2001
From: Henri Saudubray <henri.saudubray@proton.me>
Date: Mon, 12 May 2025 14:50:10 +0200
Subject: [PATCH] feat: a lot of stuff

---
 dune-project                              |  2 +-
 exm/ball.ml                               | 69 +++++++++---------
 exm/model.ml                              |  8 +++
 exm/sincos.ml                             | 30 ++++----
 exm/sqrt.ml                               |  7 +-
 exm/vdp.ml                                | 33 +++++----
 hsim.opam                                 |  2 +
 src/bin/main.ml                           | 88 ++++++++++++++---------
 src/lib/common/errors.ml                  |  2 +
 src/lib/hsim/dune                         |  2 +-
 src/lib/hsim/sim.ml                       | 75 +++++++++----------
 src/lib/hsim/solver.ml                    | 62 +++++++---------
 src/lib/hsim/statefulZ.ml                 | 34 ---------
 src/lib/hsim/types.ml                     | 60 ++++++++--------
 src/lib/solvers/csolver.ml                | 28 ++++++++
 src/lib/solvers/dune                      |  8 ++-
 src/lib/solvers/illinois.ml               | 18 ++---
 src/lib/solvers/odexx.ml                  | 17 ++---
 src/lib/{hsim => solvers}/statefulRK45.ml | 39 +++++-----
 src/lib/solvers/statefulSundials.ml       | 71 ++++++++++++++++++
 src/lib/solvers/statefulZ.ml              | 69 ++++++++++++++++++
 src/lib/solvers/zsolver.ml                | 28 ++++++++
 22 files changed, 476 insertions(+), 276 deletions(-)
 create mode 100644 exm/model.ml
 create mode 100644 src/lib/common/errors.ml
 delete mode 100644 src/lib/hsim/statefulZ.ml
 create mode 100644 src/lib/solvers/csolver.ml
 rename src/lib/{hsim => solvers}/statefulRK45.ml (65%)
 create mode 100644 src/lib/solvers/statefulSundials.ml
 create mode 100644 src/lib/solvers/statefulZ.ml
 create mode 100644 src/lib/solvers/zsolver.ml

diff --git a/dune-project b/dune-project
index 2f9f0f8..dff4a23 100644
--- a/dune-project
+++ b/dune-project
@@ -18,4 +18,4 @@
 (package
  (name hsim)
  (synopsis "An executable semantics for the simulation of hybrid systems")
- (depends ocaml))
+ (depends ocaml menhir sundialsml))
diff --git a/exm/ball.ml b/exm/ball.ml
index 9208e83..5b7eab3 100644
--- a/exm/ball.ml
+++ b/exm/ball.ml
@@ -1,6 +1,6 @@
 
 open Hsim.Types
-open Solvers
+open Solvers.Zls
 
 (* let hybrid bouncing () = (x, y) where
      rec der y  = y' init y0
@@ -11,8 +11,8 @@ open Solvers
 let of_array = Bigarray.Array1.of_array Bigarray.Float64 Bigarray.c_layout
 
 type state =
-  { zin : Zls.zarray;
-    lx : Zls.carray; (* [h';h;x';x] *)
+  { zin : zarray;
+    lx  : carray; (* [h';h;x';x] *)
     i   : bool }
 
 let g   = -9.81
@@ -20,39 +20,40 @@ let y0  = 50.0
 let y'0 = 0.0
 let x0  = 0.0
 let x'0 = 1.0
+let zsize = 1
+let csize = 4
 
-let bouncing_ball () =
-  let zfalse = Zls.zmake 1 in
-  let fder _ (y: Zls.carray) (yd: Zls.carray) =
-    if true (* y.{1} >= 0.0 *) then
-      begin yd.{0} <- g; yd.{1} <- y.{0}; yd.{2} <- 0.0; yd.{3} <- y.{2} end
-    else
-      begin yd.{0} <- 0.0; yd.{1} <- 0.0; yd.{2} <- 0.0; yd.{3} <- 0.0 end in
-  let fzer _ (y: Zls.carray) (zout: Zls.carray) = zout.{0} <- -. y.{1} in
-  let fout _ (y: Zls.carray): Zls.carray = of_array [| y.{1} |] in
-  let step ({ zin; lx; _ } as s) _ : Zls.carray * state =
-    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 } in
-  let cget s = s.lx in
-  let cset s lx = { s with lx } in
-  let zset s zin = { s with zin } in
-  let yd    = Zls.cmake 4 in
-  let zout  = Zls.cmake 1 in
-  let state = { zin=zfalse; lx=of_array [|y'0;y0;x'0;x0|]; i=true } in
-  let reset _ _ = state in
-  let jump _ = true in
-  let zsize = 1 in
-  HNode
-  { state = state;
-    fder = (fun _ _ y -> fder 0.0 y yd; yd);
-    fzer = (fun _ _ y -> fzer 0.0 y zout; zout);
-    fout = (fun s _ y -> fout s y);
-    step;
-    horizon = (fun _ -> max_float);
-    cset; cget; zset; reset; jump; zsize }
+let fder y yd =
+  if true (* y.{1} >= 0.0 *) then
+    begin yd.{0} <- g; yd.{1} <- y.{0}; yd.{2} <- 0.0; yd.{3} <- y.{2} end
+  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 jump _     = true
+let horizon _  = max_float
+let cget s     = s.lx
+let cset s lx  = { s with lx }
+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 }
+
+let bouncing_ball ()
+: (state, _, _, carray, carray, carray, zarray, carray) hnode
+= let yd         = cmake csize in
+  let zout       = cmake zsize in
+  let zfalse     = zmake 1 in
+  let fder _ _ y = fder y yd in
+  let fzer _ _ y = fzer y zout in
+  let step s _   = step s zfalse in
+  let init _     = { zin=zfalse; lx=of_array [|y'0;y0;x'0;x0|]; i=true } in
+  let reset _ _  = init () in
+  HNode { init; fder; fzer; fout; step; reset; horizon;
+          jump; cset; cget; zset; csize; zsize }
 
 let errmsg = "Too many arguments for the model (needed: 0)"
-let bouncing_ball = function
+let init = function
   | [] -> bouncing_ball ()
   | _  -> raise (Invalid_argument errmsg)
diff --git a/exm/model.ml b/exm/model.ml
new file mode 100644
index 0000000..484eeec
--- /dev/null
+++ b/exm/model.ml
@@ -0,0 +1,8 @@
+open Hsim.Types
+open Solvers.Zls
+
+module type Model =
+  sig
+    type state
+    val init : string list -> (state, 'b, 'c, carray, carray, carray, zarray, carray) hnode
+  end
diff --git a/exm/sincos.ml b/exm/sincos.ml
index 7abac36..7997d97 100644
--- a/exm/sincos.ml
+++ b/exm/sincos.ml
@@ -6,12 +6,15 @@ let of_array a = Bigarray.Array1.of_array Bigarray.Float64 Bigarray.c_layout a
 
 type state = { si : bool; sx : carray }
 
-let yd    = cmake 3
-let zout  = cmake 1
+let csize = 3
 let zsize = 1
 
-let fzer _ _ _ = zout
+let fder y yd omega =
+  yd.{0} <- omega *. y.{1}; yd.{1} <- -.omega *. y.{0}; yd.{2} <- 1.0; yd
 let fout _ _ y = of_array [| y.{0}; y.{1}; y.{2} |]
+let step { si; sx } sin0 cos0 =
+  let sx = if si then of_array [| sin0; cos0; 0.0 |] else sx in
+  of_array [| sx.{0}; sx.{1}; sx.{2} |], { sx; si=false }
 let cget s     = s.sx
 let cset s lx  = { s with sx=lx }
 let zset s _   = s
@@ -21,21 +24,20 @@ let horizon _  = max_float
 let sinus_cosinus theta0 omega =
   let sin0 = Float.sin theta0 in
   let cos0 = Float.cos theta0 in
-  let fder _ _ y =
-    yd.{0} <- omega *. y.{1}; yd.{1} <- -.omega *. y.{0}; yd.{2} <- 1.0; yd in
-  let step { si; sx } _ =
-    let sx = if si then of_array [| sin0; cos0; 0.0 |] else sx in
-    of_array [| sx.{0}; sx.{1}; sx.{2} |], { sx; si=false } in
-  let state     = { sx=of_array [| sin0; cos0; 0.0 |]; si=true } in
-  let reset _ _ = state in
-  HNode {
-    state; fder; fzer; fout; step; horizon; cset; cget; zset; zsize; reset; jump
-  }
+  let yd   = cmake csize in
+  let zout = cmake zsize in
+  let fder _ _ y = fder y yd omega in
+  let fzer _ _ _ = zout in
+  let step s _   = step s sin0 cos0 in
+  let init _     = { sx=of_array [| sin0; cos0; 0.0 |]; si=true } in
+  let reset _ _  = init () in
+  HNode { init; 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 sinus_cosinus = function
+let init = function
   | [t0; om] ->
       let t0, om = try float_of_string t0, float_of_string om
       with Failure _ -> raise (Invalid_argument errmsg_invalid) in
diff --git a/exm/sqrt.ml b/exm/sqrt.ml
index d297be8..9a86d83 100644
--- a/exm/sqrt.ml
+++ b/exm/sqrt.ml
@@ -57,6 +57,7 @@ let sqrt () =
   let yd = cmake 2 in
   let zout = cmake 1 in
   let zsize = 1 in
+  let csize = 2 in
   let s_init =
     { s_encore = false;
       s_auto = Good;
@@ -64,15 +65,15 @@ let sqrt () =
       s_zin = zmake 1 } in
   let reset _ _ = s_init in
   let jump _ = true in
-  HNode { state = s_init;
+  HNode { init = (fun _ -> s_init);
     fder = (fun s _ y -> fder s y yd; yd);
     fzer = (fun s _ y -> fzero s y zout; zout);
     fout = (fun s _ y -> fout s y);
     step = (fun s a -> fstep s a);
     horizon = (fun s -> if s.s_encore then 0.0 else max_float);
-    cset; cget; zset; zsize; reset; jump }
+    cset; cget; zset; zsize; csize; reset; jump }
 
 let errmsg = "Too many arguments to model (needed: 0)"
-let sqrt = function
+let init = function
   | [] -> sqrt ()
   | _  -> raise (Invalid_argument errmsg)
diff --git a/exm/vdp.ml b/exm/vdp.ml
index 8b22067..d224946 100644
--- a/exm/vdp.ml
+++ b/exm/vdp.ml
@@ -11,33 +11,36 @@ let mu = 5.0
 let x0 = 1.0
 let y0 = 1.0
 
+let csize = 2
 let zsize = 1
-let yd    = cmake 2
-let zout  = cmake 1
 
-let fder _ _ y =
+let fder y yd =
   yd.{0} <- y.{1};
   yd.{1} <- (mu *. (1.0 -. (y.{0} *. y.{0})) *. y.{1}) -. y.{0};
   yd
-let fzer _ _ _ = zout
 let fout _ _ y = of_array [| y.{0}; y.{1} |]
 let step { i; lx } _ =
   let lx = if i then of_array [| x0; y0 |] else lx in
   of_array [| lx.{0}; lx.{1} |], { lx; i=false }
-let cget s     = s.lx
-let cset s lx  = { s with lx }
-let zset s _   = s
-let state      = { lx=of_array [| x0; y0 |]; i=true }
-let reset _ _  = state
-let jump _     = true
-let horizon _  = max_float
+let cget s    = s.lx
+let cset s lx = { s with lx }
+let zset s _  = s
+let jump _    = true
+let horizon _ = max_float
 
-let van_der_pol () = HNode {
-  state; fder; fzer; fout; step; reset; horizon; jump; cset; cget; zset; zsize
-}
+let van_der_pol ()
+: (state, _, _, carray, carray, carray, zarray, carray) hnode
+= let yd         = cmake csize in
+  let zout       = cmake zsize in 
+  let fder _ _ y = fder y yd in
+  let fzer _ _ _ = zout in
+  let init _     = { lx=of_array [| x0; y0 |]; i=true } in
+  let reset _ _  = init () in
+  HNode { init; fder; fzer; fout; step; reset; horizon;
+          jump; cset; cget; zset; csize; zsize }
 
 let errmsg = "Too many arguments for the model (needed: 0)"
-let van_der_pol = function
+let init = function
   | [] -> van_der_pol ()
   | _  -> raise (Invalid_argument errmsg)
 
diff --git a/hsim.opam b/hsim.opam
index c68849c..bb91620 100644
--- a/hsim.opam
+++ b/hsim.opam
@@ -10,6 +10,8 @@ bug-reports: "https://codeberg.org/17maiga/hsim/issues"
 depends: [
   "dune" {>= "3.17"}
   "ocaml"
+  "menhir"
+  "sundialsml"
   "odoc" {with-doc}
 ]
 build: [
diff --git a/src/bin/main.ml b/src/bin/main.ml
index da71dba..a919e61 100644
--- a/src/bin/main.ml
+++ b/src/bin/main.ml
@@ -1,15 +1,17 @@
 
 open Hsim
-open Types
+open Solvers
 open Examples
 open Common
+open Types
 
-let sample  = ref 10
-let stop    = ref 30.0
-let greedy  = ref false
-let inplace = ref false
-let steps   = ref 1
-let model   = ref None
+let sample   = ref 10
+let stop     = ref 30.0
+let greedy   = ref false
+let inplace  = ref false
+let sundials = ref false
+let steps    = ref 1
+let model    = ref None
 
 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
@@ -21,38 +23,60 @@ let set_model s =
   | Some _ -> modelargs := s :: !modelargs
 
 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";
-  "-inplace", Arg.Set inplace,       "\tUse greedy simulation";
-  "-steps",   Arg.Int (gt0i steps),  "n \tSplit into [n] steps (default=1)";
+  "-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)";
+  "-inplace", Arg.Set inplace, "\tUse imperative solvers";
+  "-steps", Arg.Int (gt0i steps), "n \tSplit into [n] steps (default=1)";
 ]
 
 let errmsg = "Usage: " ^ Sys.executable_name ^ " [OPTIONS] MODEL\nOptions are:"
 
 let () = try Arg.parse (Arg.align opts) set_model errmsg with _ -> exit 2
 
-let output = Output.print !sample
+let m =
+  match !model with
+  | None          -> Format.eprintf "Missing model\n"; exit 2
+  | Some "ball"   -> (module Ball : Model.Model)
+  | Some "vdp"    -> (module Vdp)
+  | Some "sincos" -> (module Sincos)
+  | Some "sqrt"   -> (module Sqrt)
+  | Some s        -> Format.eprintf "Unknown model: %s\n" s; exit 2
 
-let c = StatefulRK45.(if !inplace then InPlace.csolve else Functional.csolve)
-let z = StatefulZ.(Functional.zsolve)
-let s = Solver.solver_c c z
-open Format
-let m = match !model with
-  | None          -> eprintf "Missing model\n"; exit 2
-  | Some "ball"   -> Ball.bouncing_ball
-  | Some "vdp"    -> Vdp.van_der_pol
-  | Some "sincos" -> Sincos.sinus_cosinus
-  | Some "sqrt"   -> Sqrt.sqrt
-  | Some s        -> eprintf "Unknown model: %s\n" s; exit 2
-let m = try m !modelargs with Invalid_argument s -> eprintf "%s\n" s; exit 2
+let z = if !inplace then (module StatefulZ.InPlace    : Zsolver.ZsolverC)
+                    else (module StatefulZ.Functional : Zsolver.ZsolverC)
 
-let state = if !inplace then (module State.InPlaceSimState    : State.SimState)
-                        else (module State.FunctionalSimState : State.SimState)
-let sim = if !greedy
-  then let open Sim.GreedySim(val state) in run_until_n m s
-  else let open Sim.LazySim(val state) in run_until_n m (d_of_dc s)
+let st = if !inplace then (module State.InPlaceSimState    : State.SimState)
+                     else (module State.FunctionalSimState : State.SimState)
 
-let () = sim !stop !steps output
+let () =
+  (* 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 (module InPlace    : Csolver.Csolver) *)
+  (*                         else (module Functional : Csolver.Csolver) in *)
+  (*     let open (val c) in let open (val z) in *)
+  (*     let s = Solver.solver csolve (d_of_dc zsolve) in *)
+  (*     let open Sim.LazySim(val st) in run_until_n m s *)
+  (* else *)
+  let open (val m) in
+  let m = init !modelargs in
+  let open StatefulRK45 in
+  let c = if !inplace then (module InPlace    : Csolver.CsolverC)
+          else (module Functional : Csolver.CsolverC) in
+  let open (val c) in let open (val z) in
+  let s = Solver.solver_c csolve zsolve in
+  let sim = 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) in
+  let open Solver in
+  let HNode { init; csize; zsize; fder; fzer; cget; _ } = m in
+  let state = init () in
+  let init = cget state in
+  let ivp = { size=csize; fder=(fun _ -> fder state ()); init; stop=1.0 } in
+  let zc = { size=zsize; fzer=(fun _ -> fzer state ()); init } in
+  sim !stop !steps ((), (ivp, zc)) (Output.print !sample)
 
diff --git a/src/lib/common/errors.ml b/src/lib/common/errors.ml
new file mode 100644
index 0000000..e1d582d
--- /dev/null
+++ b/src/lib/common/errors.ml
@@ -0,0 +1,2 @@
+exception TODO
+exception Internal of string
diff --git a/src/lib/hsim/dune b/src/lib/hsim/dune
index 484c26b..dd4a68b 100644
--- a/src/lib/hsim/dune
+++ b/src/lib/hsim/dune
@@ -1,3 +1,3 @@
 (library
  (name hsim)
- (libraries common solvers))
+ (libraries common))
diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml
index b0d7599..ef00fa2 100644
--- a/src/lib/hsim/sim.ml
+++ b/src/lib/hsim/sim.ml
@@ -5,13 +5,14 @@ open State
 
 module LazySim (S : SimState) =
   struct
+    module S = S
 
     (** "Lazy" 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 = S.get_init model.state solver.state in
+      (HNode model : ('ms, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
+      (DNode solver : ('ss, 'y, 'yder, 'zin, 'zout) solver)
+      : (('a, 'ms, 'ss) S.state, 'p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) lazy_sim
+    = let init (p, s) = S.get_init (model.init p) (solver.init s) in
 
       let step s i =
         let ms, ss = S.get_mstate s, S.get_sstate s in
@@ -30,10 +31,10 @@ module LazySim (S : SimState) =
                   if h <= 0.0 then S.set_mstate ms s
                   else if now >= stop then S.set_idle s
                   else if model.jump ms then
-                    let init   = model.cget ms in
+                    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 = stop -. now; init } 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;
@@ -64,32 +65,29 @@ module LazySim (S : SimState) =
         let ss = solver.reset ps (S.get_sstate s) in
         S.update ms ss (S.set_idle s) in
 
-      DNode { state; step; reset }
+      DNode { init; 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 run_on model solver input p use =
       let DNode sim = run model solver in
-      let state = match sim.step sim.state (Some input) with
-      | None, s -> s | _ -> assert false in
-      let rec loop (DNode s) =
-        let o, state = s.step s.state None in
-        match o with
-        | None   -> ()
-        | Some o -> use o; loop (DNode { s with state }) in
-      loop (DNode { sim with state })
+      let state = sim.step (sim.init p) (Some input) in
+      let state = match state with None, s -> s | _ -> assert false in
+      let rec loop state =
+        let o, state = sim.step state None in
+        match o with None -> () | Some o -> use o; loop state in
+      loop state
 
     (** Run the model on multiple inputs. *)
-    let run_on_n model solver inputs use =
-      ignore @@ List.fold_left (fun (DNode sim) i ->
-        let state = match sim.step sim.state (Some i) with
+    let run_on_n model solver inputs p use =
+      let DNode sim = run model solver in
+      ignore @@ List.fold_left (fun state i ->
+        let state = match sim.step state (Some i) with
         | None, s -> s | _ -> assert false in
-        let rec loop (DNode s) =
-          let o, state = s.step s.state None in
-          match o with
-          | None -> DNode { s with state }
-          | Some o -> use o; loop (DNode { s with state }) in
-        loop (DNode { sim with state })) (run model solver) inputs
+        let rec loop state =
+          let o, state = sim.step state None in
+          match o with None -> state | Some o -> use o; loop state in
+        loop state) (sim.init p) inputs
 
     (** Run the model autonomously until [length], or until the model stops 
         answering. *)
@@ -108,12 +106,14 @@ module LazySim (S : SimState) =
 
 module GreedySim (S : SimState) =
   struct
+    module S = 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 = S.get_init model.state solver.state in
+      (HNode model : ('ms, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
+      (DNodeC solver : ('ss, 'y, 'yder, 'zin, 'zout) solver_c)
+      : (('a, 'ms, 'ss) S.state, 'p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) greedy_sim
+    = let init (m, s) = S.get_init (model.init m) (solver.init s) in
 
       let rec step s i =
         let ms, ss = S.get_mstate s, S.get_sstate s in
@@ -132,7 +132,7 @@ module GreedySim (S : SimState) =
                 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 } 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; 
@@ -171,20 +171,21 @@ module GreedySim (S : SimState) =
         let ss = solver.reset sp (S.get_sstate s) in
         S.update ms ss (S.set_idle s) in
 
-      DNode { state; step; reset }
+      DNode { init; 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 run_on model solver input p use =
       let DNode sim = run model solver in
-      let o, _ = sim.step sim.state input in
+      let o, _ = sim.step (sim.init p) input in
       List.iter use o
 
     (** Run the model on multiple inputs. *)
-    let run_on_n model solver inputs use =
-      let o, _ = List.fold_left (fun (acc, DNode sim) i ->
-        let o, state = sim.step sim.state i in
-        o::acc, DNode { sim with state }) ([], run model solver) inputs in
+    let run_on_n model solver inputs p 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.init p) inputs in
       List.iter use (List.concat (List.rev o))
 
     (** Run the model autonomously until [length], or until the model stops 
diff --git a/src/lib/hsim/solver.ml b/src/lib/hsim/solver.ml
index c8be4a7..e427b46 100644
--- a/src/lib/hsim/solver.ml
+++ b/src/lib/hsim/solver.ml
@@ -5,7 +5,8 @@ open Types
 type ('y, 'yder) ivp =
   { init : 'y;                         (** [y₀]: initial value of y.          *)
     fder : time -> 'y -> 'yder;        (** [dy/dt]: derivative of y.          *)
-    stop : time }                      (** Stop time.                         *)
+    stop : time;                       (** Stop time.                         *)
+    size : int }
 
 (** A zero-crossing expression. *)
 type ('y, 'zout) zc =
@@ -17,72 +18,62 @@ type ('y, 'zout) zc =
   - an initial value problem as parameter;
   - an horizon to reach as input;
   - an actual time reached and dense solution as output *)
-type ('y, 'yder) csolver =
-  (('y, 'yder) ivp, time, time * (time -> 'y)) dnode
+type ('s, 'y, 'yder) csolver =
+  ('s, ('y, 'yder) ivp, time, time * (time -> 'y)) dnode
 
 (** An ODE solver can optionally provide a state copy method, in which case
     greedy simulation is possible. *)
-type ('y, 'yder) csolver_c =
-  (('y, 'yder) ivp, time, time * (time -> 'y)) dnode_c
+type ('s, 'y, 'yder) csolver_c =
+  ('s, ('y, 'yder) ivp, time, time * (time -> 'y)) dnode_c
 
 (** A zero-crossing solver is a synchronous function with:
   - a zero-crossing expression as parameter;
   - a time and dense solution as input;
   - an actual time reached and optional zero-crossing as output *)
-type ('y, 'zin, 'zout) zsolver =
-  (('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode
+type ('s, 'y, 'zin, 'zout) zsolver =
+  ('s, ('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode
 
 (** A zero-crossing solver can optionally provide a state copy method, in which
     case greedy simulation is possible. *)
-type ('y, 'zin, 'zout) zsolver_c =
-  (('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode_c
+type ('s, 'y, 'zin, 'zout) zsolver_c =
+  ('s, ('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode_c
 
 (** A solver is a synchronous function with:
   - an initial value problem and zero-crossing expression as parameter;
   - an horizon to reach as input;
   - an actual time, dense solution and optional zero-crossing as output *)
-type ('y, 'yder, 'zin, 'zout) solver =
-  (('y, 'yder) ivp * ('y, 'zout) zc,
+type ('s, 'y, 'yder, 'zin, 'zout) solver =
+  ('s,
+   ('y, 'yder) ivp * ('y, 'zout) zc,
    time,
    time * (time -> 'y) * 'zin option) dnode
 
 (** A solver can optionally provide a state copy method, in which case greedy
     simulation is possible. *)
-type ('y, 'yder, 'zin, 'zout) solver_c =
-  (('y, 'yder) ivp * ('y, 'zout) zc,
+type ('s, 'y, 'yder, 'zin, 'zout) solver_c =
+  ('s,
+   ('y, 'yder) ivp * ('y, 'zout) zc,
    time,
    time * (time -> 'y) * 'zin option) dnode_c
 
-let csolver_from_c (DNodeC csolver : ('y, 'yder) csolver_c)
-  : ('y, 'yder) csolver
-= DNode { state = csolver.state; step = csolver.step; reset = csolver.reset }
-
-let zsolver_from_c (DNodeC zsolver : ('y, 'zin, 'zout) zsolver_c)
-  : ('y, 'zin, 'zout) zsolver
-= DNode { state = zsolver.state; step = zsolver.step; reset = zsolver.reset }
-
-let solver_from_c (DNodeC solver : ('y, 'yder, 'zin, 'zout) solver_c)
-  : ('y, 'yder, 'zin, 'zout) solver
-= DNode { state = solver.state; step = solver.step; reset = solver.reset }
-
 (** Build a full solver from an ODE solver and a zero-crossing solver. *)
-let solver (DNode csolver : ('y, 'yder) csolver)
-           (DNode zsolver : ('y, 'zin, 'zout) zsolver)
-    : ('y, 'yder, 'zin, 'zout) solver =
-  let state = csolver.state, zsolver.state in
+let solver (DNode csolver : ('sc, 'y, 'yder) csolver)
+           (DNode zsolver : ('sz, 'y, 'zin, 'zout) zsolver)
+    : ('sc * 'sz, 'y, 'yder, 'zin, 'zout) solver =
+  let init (ivp, zc) = csolver.init ivp, zsolver.init zc in
   let step (cstate, zstate) h =
     let (h, f), cstate = csolver.step cstate h in
     let (h, z), zstate = zsolver.step zstate (h, f) in
     (h, f, z), (cstate, zstate) in
   let reset (ivp, zc) (cstate, zstate) =
     csolver.reset ivp cstate, zsolver.reset zc zstate in
-  DNode { state; step; reset }
+  DNode { init ; step; reset }
 
 (** Build a full solver supporting state copies. *)
-let solver_c (DNodeC csolver : ('y, 'yder) csolver_c)
-             (DNodeC zsolver : ('y, 'zin, 'zout) zsolver_c)
-    : ('y, 'yder, 'zin, 'zout) solver_c =
-  let state = csolver.state, zsolver.state in
+let solver_c (DNodeC csolver : ('sc, 'y, 'yder) csolver_c)
+             (DNodeC zsolver : ('sz, 'y, 'zin, 'zout) zsolver_c)
+    : ('sc * 'sz, 'y, 'yder, 'zin, 'zout) solver_c =
+  let init (ivp, zc) = csolver.init ivp, zsolver.init zc in
   let step (cstate, zstate) h =
     let (h, f), cstate = csolver.step cstate h in
     let (h, z), zstate = zsolver.step zstate (h, f) in
@@ -91,5 +82,4 @@ let solver_c (DNodeC csolver : ('y, 'yder) csolver_c)
     csolver.reset ivp cstate, zsolver.reset zc zstate in
   let copy (cstate, zstate) =
     csolver.copy cstate, zsolver.copy zstate in
-  DNodeC { state; step; reset; copy }
-
+  DNodeC { init; step; reset; copy }
diff --git a/src/lib/hsim/statefulZ.ml b/src/lib/hsim/statefulZ.ml
deleted file mode 100644
index cf8787c..0000000
--- a/src/lib/hsim/statefulZ.ml
+++ /dev/null
@@ -1,34 +0,0 @@
-
-open Types
-open Solvers
-open Solver
-
-module Functional =
-  struct
-
-    type ('state, 'vec) state = { state: 'state; vec: 'vec }
-
-    let zsolve : (Zls.carray, Zls.zarray, Zls.carray) zsolver_c =
-      let state =
-        { state = Illinois.initialize 0 (fun _ _ _ -> ()) (Zls.cmake 0);
-          vec = Zls.zmake 0 } in
-      let reset { fzer; init; size } { vec; _ } =
-        let fzer t cvec zout = let zout' = fzer t cvec in Zls.blit zout' zout in
-        { state = Illinois.initialize size fzer init;
-          vec = if Zls.length vec = size then vec else Zls.zmake size } in
-
-      let step ({ state; vec } as s) (h, fder) =
-        let y1 = fder h in
-        let fder h _ = let y = fder h in Zls.blit y y1 in
-        Illinois.step state h y1;
-        let v = Illinois.has_roots state in
-        if v then
-          let h = Illinois.find state (fder, y1) vec in
-          (h, Some vec), s
-        else (h, None), s in
-
-      let copy s = s in
-
-      DNodeC { state; step; reset; copy }
-
-  end
diff --git a/src/lib/hsim/types.ml b/src/lib/hsim/types.ml
index e8fdaa7..0054331 100644
--- a/src/lib/hsim/types.ml
+++ b/src/lib/hsim/types.ml
@@ -14,21 +14,21 @@ type 'a value =
 type 'a signal = 'a value option
 
 (** A discrete node. *)
-type ('p, 'a, 'b) dnode =
+type ('s, 'p, 'a, 'b) dnode =
   DNode :
-    { state : 'ds;
-      step  : 'ds -> 'a -> 'b * 'ds;
-      reset : 'p -> 'ds -> 'ds;
-    } -> ('p, 'a, 'b) dnode
+    { init  : 'p -> 's;
+      step  : 's -> 'a -> 'b * 's;
+      reset : 'p -> 's -> 's;
+    } -> ('s, 'p, 'a, 'b) dnode
 
 (** A discrete node which supports a state copy. *)
-type ('p, 'a, 'b) dnode_c =
+type ('s, 'p, 'a, 'b) dnode_c =
   DNodeC :
-    { state : 'ds;
-      step  : 'ds -> 'a -> 'b * 'ds;
-      reset : 'p -> 'ds -> 'ds;
-      copy  : 'ds -> 'ds;
-    } -> ('p, 'a, 'b) dnode_c
+    { init  : 'p -> 's;
+      step  : 's -> 'a -> 'b * 's;
+      reset : 'p -> 's -> 's;
+      copy  : 's -> 's;
+    } -> ('s, 'p, 'a, 'b) dnode_c
 
 (** A continuous node. *)
 type ('a, 'b, 'y, 'yder) cnode =
@@ -39,33 +39,33 @@ type ('a, 'b, 'y, 'yder) cnode =
     } -> ('a, 'b, 'y, 'yder) cnode
 
 (** A hybrid node. *)
-type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
+type ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
   HNode :
-    { state : 'hs;
-      step : 'hs -> 'a -> 'b * 'hs;    (** Discrete step function.            *)
-      fder : 'hs -> 'a -> 'y -> 'yder; (** Continuous derivative function.    *)
-      fout : 'hs -> 'a -> 'y -> 'b;    (** Continuous output function.        *)
-      fzer : 'hs -> 'a -> 'y -> 'zout; (** Continuous zero-crossing function. *)
-      reset : 'p -> 'hs -> 'hs;        (** Reset function.                    *)
-      horizon : 'hs -> time;           (** Next integration horizon.          *)
-      jump : 'hs -> bool;              (** Discontinuity flag.                *)
-      cget : 'hs -> 'y;                (** Get continuous state.              *)
-      cset : 'hs -> 'y -> 'hs;         (** Set continuous state.              *)
-      zset : 'hs -> 'zin -> 'hs;       (** Set zero-crossing state.           *)
+    { init : 'p -> '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
+    } -> ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
 
 (** 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 ('s, 'p, 'a, 'b) lazy_sim =
+  ('s, '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
+type ('s, 'p, 'a, 'b) greedy_sim =
+  ('s, 'p, 'a value, 'b value list) dnode
 
 (** Utils *)
 
-let d_of_dc (DNodeC { state; step; reset; _ }) =
-  DNode { state; step; reset }
+let d_of_dc (DNodeC { init; step; reset; _ }) = DNode { init; step; reset }
diff --git a/src/lib/solvers/csolver.ml b/src/lib/solvers/csolver.ml
new file mode 100644
index 0000000..71002b1
--- /dev/null
+++ b/src/lib/solvers/csolver.ml
@@ -0,0 +1,28 @@
+
+open Hsim.Types
+open Hsim.Solver
+open Zls
+
+module type Csolver =
+  sig
+    type ('a, 'b) state
+    type session
+    type vec
+    val csolve : ((session, vec) state, carray, carray) csolver
+  end
+
+module type CsolverC =
+  sig
+    type ('a, 'b) state
+    type session
+    type vec
+    val csolve : ((session, vec) state, carray, carray) csolver_c
+  end
+
+module CsolverOfC =
+  functor (S : CsolverC) -> (struct
+    type ('a, 'b) state = ('a, 'b) S.state
+    type session = S.session
+    type vec = S.vec
+    let csolve = d_of_dc S.csolve
+  end : Csolver)
diff --git a/src/lib/solvers/dune b/src/lib/solvers/dune
index b92f525..099ae6f 100644
--- a/src/lib/solvers/dune
+++ b/src/lib/solvers/dune
@@ -1,4 +1,8 @@
-(env (dev (flags (:standard -w -9-27-32))))
+(env (dev (flags (:standard -w -32))))
 
 (library
-  (name solvers))
+  (name solvers)
+  (libraries
+    hsim
+    ;sundialsml
+    ))
diff --git a/src/lib/solvers/illinois.ml b/src/lib/solvers/illinois.ml
index 2919029..434be45 100644
--- a/src/lib/solvers/illinois.ml
+++ b/src/lib/solvers/illinois.ml
@@ -27,7 +27,7 @@ let get_check_root rdir =
   let check_down   x0 x1 = if x0 > 0.0 && x1 <= 0.0 then -1l else 0l in
   let check_either x0 x1 = if x0 < 0.0 && x1 >= 0.0 then  1l else
                            if x0 > 0.0 && x1 <= 0.0 then -1l else 0l in
-  let no_check x0 x1 = 0l in
+  let no_check _x0 _x1 = 0l in
 
   match rdir with
   | Up     -> check_up
@@ -118,7 +118,7 @@ type t = {
   }
 
 (* Called from find when bothf_valid = false to initialise f1. *)
-let reinitialize ({ g; f1 = f1; t1 = t1 } as s) t c =
+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;
@@ -148,10 +148,10 @@ let initialize nroots g c =
   s
 
 
-let num_roots { f0 } = Zls.length f0
+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 =
+let step ({ g; f0 = f0; f1 = f1; t1 = t1; _ } as s) t c =
   (* swap f0 and f1; f0 takes the previous value of f1 *)
   s.f0 <- f1;
   s.t0 <- t1;
@@ -184,7 +184,7 @@ let resolve_intervals r1 r2 =
                    (possible) zero-crossing in (f_mid, f_right]
  *)
 let check_interval calc_zc f_left f_mid =
-  let check i r x0 x1 =
+  let check _i r x0 x1 =
     let rv = calc_zc x0 x1 in
     let r' = if rv = 0l then SearchRight
              else if x1 = 0.0 then FoundMid
@@ -340,17 +340,17 @@ let find s (dky, c) roots = find s (dky, c) roots
 
 (* is there a root? [has_root s: bool] is true is there is a change in sign *)
 (* for one component [i in [0..length f0 - 1]] beetwen [f0.(i)] and [f1.(i)] *)
-let has_roots { bothf_valid = bothf_valid; t0; f0; t1; f1; calc_zc = calc_zc }
-  = bothf_valid && (check_interval calc_zc f0 f1 <> SearchRight)
+let has_roots { bothf_valid; f0; f1; calc_zc; _ } =
+  bothf_valid && (check_interval calc_zc f0 f1 <> SearchRight)
 
-let takeoff { bothf_valid = bothf_valid; f0; f1 } =
+let takeoff { bothf_valid; f0; f1; _ } =
   bothf_valid && (takeoff f0 f1)
 
 (* returns true if a signal has moved from zero to a stritly positive value *)
 (* Added by MP. Ask Tim if this code is necessary, that is, what happens *)
 (* with function [find] when the signal is taking off from [0.0] to a *)
 (* strictly positive value *)
-let find_takeoff ({ f0; f1 } as s) roots =
+let find_takeoff ({ f0; f1; _ } as s) roots =
   let calc_zc x0 x1 =
     if (x0 = 0.0) && (x1 > 0.0) then 1l else 0l in
   let b = update_roots calc_zc f0 f1 roots in
diff --git a/src/lib/solvers/odexx.ml b/src/lib/solvers/odexx.ml
index 3eca888..31637aa 100644
--- a/src/lib/solvers/odexx.ml
+++ b/src/lib/solvers/odexx.ml
@@ -156,7 +156,7 @@ struct (* {{{1 *)
   (* NB: y must be the initial state vector (y_0)
    *     k(0) must be the initial deriviatives vector (dy_0) *)
   let initial_stepsize { initial_step_size; abs_tol; rel_tol; max_step;
-                         time; y; hmax; k } =
+                         time; y; hmax; k; _ } =
     let hmin = 16.0 *. epsilon_float *. abs_float time in
     match initial_step_size with
     | Some h -> minmax hmin max_step h
@@ -168,7 +168,8 @@ struct (* {{{1 *)
       in
       max hmin (if hmax *. rh > 1.0 then 1.0 /. rh else hmax)
 
-  let reinitialize ?rhsfn ({ stop_time; min_step; max_step; sysf } as s) t ny =
+  let reinitialize
+      ?rhsfn ({ stop_time; min_step; max_step; sysf; _ } as s) t ny =
     Bigarray.Array1.blit ny s.y;
     s.time      <- t;
     s.last_time <- t;
@@ -250,9 +251,9 @@ struct (* {{{1 *)
 
   (* TODO: add stats: nfevals, nfailed, nsteps *)
   let step s t_limit user_y =
-    let { stop_time; min_step; abs_tol; rel_tol;
+    let { stop_time; abs_tol; rel_tol;
           sysf = f; time = t; h = h; hmax = hmax;
-          k = k; y = y; yold = ynew; } = s in
+          k = k; y = y; yold = ynew; _ } = s in
 
     (* First Same As Last (FSAL) swap; doing it after the previous
        step invalidates the interpolation routine. *)
@@ -323,7 +324,7 @@ struct (* {{{1 *)
     s.h <- nexth;
     s.time
 
-  let get_dky { last_time = t; time = t'; h = h; yold = y; k = k } yi ti kd =
+  let get_dky { last_time = t; time = t'; yold = y; k; _ } yi ti kd =
 
     if kd > 0 then
       failwith
@@ -355,11 +356,11 @@ struct (* {{{1 *)
     done
 
   (* copy functions *)
-  let copy ({ last_time; time; h; yold; k } as s) =
+  let copy ({ last_time; time; h; yold; k; _ } as s) =
     { s with last_time; time; h; yold = Zls.copy yold; k = Zls.copy_matrix k }
 
-  let blit { last_time = l1; time = t1; h = h1; yold = yhold1; k = k1 }
-      ({ last_time; time; h; yold; k } as s2) =
+  let blit { last_time = l1; time = t1; yold = yhold1; k = k1; _ }
+      ({ yold; k; _ } as s2) =
     s2.last_time <- l1; s2.time <- t1;
     Zls.blit yhold1 yold; Zls.blit_matrix k1 k
 
diff --git a/src/lib/hsim/statefulRK45.ml b/src/lib/solvers/statefulRK45.ml
similarity index 65%
rename from src/lib/hsim/statefulRK45.ml
rename to src/lib/solvers/statefulRK45.ml
index 2689681..6a050eb 100644
--- a/src/lib/hsim/statefulRK45.ml
+++ b/src/lib/solvers/statefulRK45.ml
@@ -1,25 +1,24 @@
 
-open Types
-open Solvers
-open Solver
+open Hsim.Types
+open Hsim.Solver
+open Zls
 
-module Functional =
+module Functional : Csolver.CsolverC =
   struct
     type ('state, 'vec) state = { state: 'state; vec: 'vec }
+    type session = Odexx.Ode45.t
+    type vec = carray
 
-    let csolve : (Zls.carray, Zls.carray) csolver_c =
+    let csolve : ((session, vec) state, carray, carray) csolver_c =
       let open Odexx.Ode45 in
 
-      let state =
+      let init _ =
         let v = Zls.cmake 0 in
         let state = initialize (fun _ _ _ -> ()) (vec v) in
         set_stop_time state 1.0; { state; vec=v } in
 
-      let reset
-        ({ fder; init; stop }: (Zls.carray, Zls.carray) ivp)
-        (_: (t, Zls.carray) state)
-        : (t, Zls.carray) state
-      = let fder t cvec dvec = Zls.blit (fder t cvec) dvec in
+      let reset { fder; init; stop; _ } _ = 
+        let fder t cvec dvec = Zls.blit (fder t cvec) dvec in
         let state = initialize fder (vec init) in
         set_stop_time state stop;
         { state; vec = init } in
@@ -33,25 +32,25 @@ module Functional =
 
       let copy { state; vec } = { state; vec } in
 
-      DNodeC { state; step; reset; copy }
+      DNodeC { init; step; reset; copy }
   end
 
-module InPlace =
+module InPlace : Csolver.CsolverC =
   struct
+    type ('state, 'vec) state = { mutable state: 'state; mutable vec : 'vec }
+    type session = Odexx.Ode45.t
+    type vec = carray
 
-    type ('state, 'vec) state =
-      { mutable state: 'state; mutable vec : 'vec }
-
-    let csolve : (Zls.carray, Zls.carray) csolver_c =
+    let csolve : ((session, vec) state, carray, carray) csolver_c =
       let open Odexx.Ode45 in
 
-      let state =
+      let init _ =
         let v = Zls.cmake 0 in
         let state = initialize (fun _ _ _ -> ()) (vec v) in
         set_stop_time state 1.0;
         { state; vec=v } in
 
-      let reset { fder: time -> Zls.carray -> Zls.carray; init; stop } s =
+      let reset { fder; init; stop; _ } s =
         let fder t cvec dvec =
           let dvec' = fder t cvec in Zls.blit dvec' dvec in
         let state = initialize fder (vec init) in
@@ -66,5 +65,5 @@ module InPlace =
       let copy { state; vec } =
         { state = copy state; vec = Zls.copy vec } in
 
-      DNodeC { state; reset; step; copy }
+      DNodeC { init; reset; step; copy }
   end
diff --git a/src/lib/solvers/statefulSundials.ml b/src/lib/solvers/statefulSundials.ml
new file mode 100644
index 0000000..5308f7e
--- /dev/null
+++ b/src/lib/solvers/statefulSundials.ml
@@ -0,0 +1,71 @@
+(*
+open Hsim.Types
+open Hsim.Solver
+open Zls
+
+module Functional : Csolver.Csolver =
+  struct
+    type ('state, 'vec) state = { state : 'state; vec : 'vec }
+    type session = (Sundials_RealArray.t, Nvector_serial.kind) Cvode.session
+    type vec = carray
+
+    let csolve : ((session, vec) state, carray, carray) csolver =
+      let open Cvode in
+
+      let init { size; fder=_; _ } =
+        let vec = cmake size in
+        let state = init Adams default_tolerances (fun _ _ _ -> ()) 0.
+                         (Nvector_serial.wrap vec) in
+        set_stop_time state 1.0;
+        { state; vec } in
+
+      let reset { init=i; fder; stop; _ } { vec; _ } =
+        let fder t cvec dvec =
+          let dvec' = fder t cvec in blit dvec' dvec in
+        let state =
+          Cvode.init Adams default_tolerances fder 0. (Nvector_serial.wrap i) in
+        set_stop_time state stop;
+        { state; vec } in
+
+      let step ({ state; vec } as s) h =
+        let y = Nvector_serial.wrap vec in
+        let h, _ = solve_one_step state h y in
+        let f t = get_dky state y t 0; Nvector_serial.unwrap y in
+        (h, f), s in
+
+      DNode { init; reset; step }
+  end
+
+module InPlace : Csolver.Csolver =
+  struct
+    type ('state, 'vec) state = { mutable state: 'state; mutable vec : 'vec }
+
+    type session = (Sundials_RealArray.t, Nvector_serial.kind) Cvode.session
+    type vec = carray
+
+    let csolve : ((session, vec) state, carray, carray) csolver =
+      let open Cvode in
+
+      let init { size; fder=_; _ } =
+        let vec = cmake size in
+        let state = init Adams default_tolerances (fun _ _ _ -> ()) 0.
+                         (Nvector_serial.wrap vec) in
+        set_stop_time state 1.0;
+        { state; vec } in
+
+      let reset { init=i; fder; _ } s =
+        let fder t cvec dvec =
+          let dvec' = fder t cvec in blit dvec' dvec in
+        let state =
+          Cvode.init Adams default_tolerances fder 0. (Nvector_serial.wrap i) in
+        set_stop_time state 1.0; s.state <- state; s.vec <- i; s in
+
+      let step s h =
+        let y = Nvector_serial.wrap s.vec in
+        let h, _ = solve_one_step s.state h y in
+        let f t = get_dky s.state y t 0; Nvector_serial.unwrap y in
+        (h, f), s in 
+
+      DNode { init; reset; step }
+  end
+*)
diff --git a/src/lib/solvers/statefulZ.ml b/src/lib/solvers/statefulZ.ml
new file mode 100644
index 0000000..bfb83a1
--- /dev/null
+++ b/src/lib/solvers/statefulZ.ml
@@ -0,0 +1,69 @@
+
+open Hsim.Types
+open Hsim.Solver
+open Zls
+
+module Functional : Zsolver.ZsolverC =
+  struct
+
+    type ('state, 'vec) state = { state: 'state; vec: 'vec }
+    type session = Illinois.t
+    type vec = zarray
+
+    let zsolve : ((session, vec) state, carray, vec, carray) zsolver_c =
+      let open Illinois in
+
+      let init _ =
+        { state = initialize 0 (fun _ _ _ -> ()) (cmake 0);
+          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
+        { state = initialize size fzer init;
+          vec = if length vec = size then vec else zmake size } in
+
+      let step ({ state; vec } as s) (h, fder) =
+        let y1 = fder h in
+        let fder h _ = let y = fder h in blit y y1 in
+        step state h y1;
+        if has_roots state then
+          let h = find state (fder, y1) vec in
+          (h, Some vec), s
+        else (h, None), s in
+
+      let copy s = s in
+
+      DNodeC { init; step; reset; copy }
+  end
+
+module InPlace : Zsolver.ZsolverC =
+  struct
+    type ('state, 'vec) state = { mutable state : 'state; mutable vec : 'vec }
+    type session = Illinois.t
+    type vec = zarray
+
+    let zsolve : ((session, vec) state, carray, vec, carray) zsolver_c =
+      let open Illinois in
+
+      let init _ =
+        { state=initialize 0 (fun _ _ _ -> ()) (cmake 0);
+          vec=zmake 0 } in
+
+      let reset { size; init; fzer } s =
+        let fzer t cvec zout = let zout' = fzer t cvec in blit zout' zout in
+        s.state <- initialize size fzer init;
+        if length s.vec <> size then s.vec <- zmake size; s in
+
+      let step ({ state; vec } as s) (h, fder) =
+        let y = fder h in
+        let fder h _ = let y' = fder h in blit y' y in
+        step state h y;
+        if has_roots state then
+          let h = find state (fder, y) vec in
+          (h, Some vec), s
+        else (h, None), s in
+
+      let copy _ = raise Common.Errors.TODO in
+      
+      DNodeC { init; step; reset; copy }
+  end
diff --git a/src/lib/solvers/zsolver.ml b/src/lib/solvers/zsolver.ml
new file mode 100644
index 0000000..2d63d4b
--- /dev/null
+++ b/src/lib/solvers/zsolver.ml
@@ -0,0 +1,28 @@
+
+open Hsim.Types
+open Hsim.Solver
+open Zls
+
+module type Zsolver =
+  sig
+    type ('a, 'b) state
+    type session
+    type vec
+    val zsolve : ((session, vec) state, carray, zarray, carray) zsolver
+  end
+
+module type ZsolverC =
+  sig
+    type ('a, 'b) state
+    type session
+    type vec
+    val zsolve : ((session, vec) state, carray, zarray, carray) zsolver_c
+  end
+
+module ZsolverOfC =
+  functor (S : ZsolverC) -> (struct
+    type ('a, 'b) state = ('a, 'b) S.state
+    type session = S.session
+    type vec = S.vec
+    let zsolve = d_of_dc S.zsolve
+  end : Zsolver)