diff --git a/.ocamlformat b/.ocamlformat
new file mode 100644
index 0000000..5143690
--- /dev/null
+++ b/.ocamlformat
@@ -0,0 +1,2 @@
+version=0.27.0
+margin=80
diff --git a/bin/main.ml b/bin/main.ml
deleted file mode 100644
index 7bf6048..0000000
--- a/bin/main.ml
+++ /dev/null
@@ -1 +0,0 @@
-let () = print_endline "Hello, World!"
diff --git a/lib/dune b/lib/dune
deleted file mode 100644
index b00369f..0000000
--- a/lib/dune
+++ /dev/null
@@ -1,2 +0,0 @@
-(library
- (name hsim))
diff --git a/bin/dune b/src/bin/dune
similarity index 100%
rename from bin/dune
rename to src/bin/dune
diff --git a/src/bin/main.ml b/src/bin/main.ml
new file mode 100644
index 0000000..15659b5
--- /dev/null
+++ b/src/bin/main.ml
@@ -0,0 +1,4 @@
+open Hsim.Types
+
+let _x : 'a value = { start = 0.; length = 0.; u = (fun _ -> 0) }
+let () = print_endline "Hello, World!"
diff --git a/src/lib/common/dune b/src/lib/common/dune
new file mode 100644
index 0000000..4ce44f9
--- /dev/null
+++ b/src/lib/common/dune
@@ -0,0 +1,2 @@
+(library
+  (name common))
diff --git a/src/lib/common/monad.ml b/src/lib/common/monad.ml
new file mode 100644
index 0000000..e614a29
--- /dev/null
+++ b/src/lib/common/monad.ml
@@ -0,0 +1,26 @@
+
+module type Monad = sig
+  type 'a t
+
+  val return : 'a -> 'a t
+  val bind   : 'a t -> ('a -> 'b t) -> 'b t
+end
+
+module type FullMonad = sig
+  type 'a t
+
+  val return : 'a -> 'a t
+  val bind   : 'a t -> ('a -> 'b t) -> 'b t
+  val (>>=)  : 'a t -> ('a -> 'b t) -> 'b t
+  val (let*) : 'a t -> ('a -> 'b t) -> 'b t
+  val join   : 'a t t -> 'a t
+end
+
+module Expand (M : Monad) = struct
+  include M
+
+  let (>>=)  = M.bind
+  let (let*) = M.bind
+  let join m = M.bind m (fun m -> m)
+end
+
diff --git a/src/lib/common/monad.mli b/src/lib/common/monad.mli
new file mode 100644
index 0000000..2543fa3
--- /dev/null
+++ b/src/lib/common/monad.mli
@@ -0,0 +1,20 @@
+
+module type Monad = sig
+  type 'a t
+
+  val return : 'a -> 'a t
+  val bind   : 'a t -> ('a -> 'b t) -> 'b t
+end
+
+module type FullMonad = sig
+  type 'a t
+
+  val return : 'a -> 'a t
+  val bind   : 'a t -> ('a -> 'b t) -> 'b t
+  val (>>=)  : 'a t -> ('a -> 'b t) -> 'b t
+  val (let*) : 'a t -> ('a -> 'b t) -> 'b t
+  val join   : 'a t t -> 'a t
+end
+
+module Expand : functor (M : Monad) -> FullMonad with type 'a t = 'a M.t
+
diff --git a/src/lib/common/state.ml b/src/lib/common/state.ml
new file mode 100644
index 0000000..a09b07a
--- /dev/null
+++ b/src/lib/common/state.ml
@@ -0,0 +1,23 @@
+
+open Utils
+
+module Make (S : sig
+  type t
+end) =
+  struct
+
+    module State = struct
+      type 'a t = S.t -> 'a * S.t
+
+      let return   = pair
+      let bind m f = uncurry f @. m
+    end
+
+    module M = Monad.Expand (State)
+    include M
+
+    let get () s = s, s
+    let set x _  = (), x
+    let run m    = fst @. m
+  end
+
diff --git a/src/lib/common/state.mli b/src/lib/common/state.mli
new file mode 100644
index 0000000..dfb6289
--- /dev/null
+++ b/src/lib/common/state.mli
@@ -0,0 +1,17 @@
+
+module Make (S : sig
+  type t
+end) : sig
+  type 'a t
+
+  val return : 'a -> 'a t
+  val bind   : 'a t -> ('a -> 'b t) -> 'b t
+  val (>>=)  : 'a t -> ('a -> 'b t) -> 'b t
+  val (let*) : 'a t -> ('a -> 'b t) -> 'b t
+
+  val get    : unit -> S.t t
+  val set    : S.t -> unit t
+
+  val run    : 'a t -> S.t -> 'a
+end
+
diff --git a/src/lib/common/utils.ml b/src/lib/common/utils.ml
new file mode 100644
index 0000000..bf1126f
--- /dev/null
+++ b/src/lib/common/utils.ml
@@ -0,0 +1,7 @@
+
+let pair = fun a b -> a, b
+
+let uncurry = fun f (a, b) -> f a b
+
+let (@.) = fun f g x -> f @@ g x
+
diff --git a/src/lib/hsim/dune b/src/lib/hsim/dune
new file mode 100644
index 0000000..b810b4e
--- /dev/null
+++ b/src/lib/hsim/dune
@@ -0,0 +1,4 @@
+(library
+ (name hsim)
+ (libraries common)
+ (modules_without_implementation types zls))
diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml
new file mode 100644
index 0000000..fcf2d53
--- /dev/null
+++ b/src/lib/hsim/sim.ml
@@ -0,0 +1,98 @@
+
+open Types
+
+(** Step mode: discrete or continuous *)
+type mode =
+  | Discrete
+  | Continuous
+
+(** Simulation status:
+  - [Idle]:    Waiting for input, no activity;
+  - [Running]: Currently integrating: step [mode], current [input], [now] 
+               timestamp, and [stop] time. *)
+type 'a status =
+  | Idle    : 'a status
+  | Running :
+      { mode  : mode;                        (** Step mode.                   *)
+        input : 'a value;                    (** Function to integrate.       *)
+        now   : time;                        (** Current time of integration. *)
+        stop  : time;                        (** How long until we stop.      *)
+      } -> 'a status
+
+(** Internal state of the simulation node: model state, solver state and current
+    simulation status. *)
+type ('a, 'ms, 'ss) state =
+  { status : 'a status;                      (** Current simulation status.   *)
+    mstate : 'ms;                            (** Model state.                 *)
+    sstate : 'ss }                           (** Solver state.                *)
+
+(** Offset the [input] function by [now]. *)
+let offset (input: 'a value) (now: time) : time -> 'a =
+  fun t -> input.u ((now -. input.start) +. t)
+
+let simulate (HNode model: ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
+             (DNode solver : ('y, 'yder, 'zin, 'zout) solver)
+                 : ('p, 'a signal, 'b signal) dnode =
+  let s = { status = Idle; mstate = model.s; sstate = solver.s } in 
+  let step (state: ('a, _, _) state) (input: 'a signal)
+      : 'b signal * ('a, _, _) state =
+    match input, state.status with
+    | Some i, _ ->
+        (* New input. *)
+        let status = Running
+          { mode = Continuous; input = i; now = 0.0; stop = i.length } in
+        None, { state with status }
+    | None, Idle ->
+        (* Waiting for input. *)
+        None, state
+    | None, Running ({ mode = Discrete; _ } as r) ->
+        (* Discrete step. *)
+        let o, mstate = model.step state.mstate (r.input.u r.now) in
+        (* Four possible outcomes: *)
+        let state =
+          let h = model.horizon mstate in
+          if h <= 0.0 then                   (* - Cascade (new discrete step) *)
+            { state with mstate }
+          else if r.now >= r.stop then       (* - Reached horizon             *)
+            s
+          else if model.jump mstate then     (* - Discontinuity: reset solver *)
+            let y      = model.cget mstate in
+            let fder t = model.fder mstate @@ offset r.input r.now t in
+            let fzer t = model.fzer mstate @@ offset r.input r.now t in
+            let ivp    = { fder; stop = r.stop -. r.now; init = y } in
+            let zc     = { yc = y; fzer } in
+            let sstate = solver.reset (ivp, zc) state.sstate in
+            let status = Running { r with mode = Continuous } in
+            { status; mstate; sstate }
+          else                               (* - Continue                    *)
+            { state with status = Running { r with mode = Continuous } }
+        in Some { start = r.now; length = 0.0; u = fun _ -> o }, state
+    | None, Running ({ mode = Continuous; _ } as r) ->
+        (* Continuous step *)
+        let (h, f, z), sstate = solver.step state.sstate r.stop in
+        let mstate = model.cset state.mstate (f h) in
+        (* Three possible outcomes: *)
+        let now = r.input.start +. h in
+        let state = match z with
+        | None ->
+            let status =
+              if h >= r.stop then                 (* Reached the end.         *)
+                Running { r with mode = Discrete; now }
+              else                                (* Not finished integrating *)
+                Running { r with now } in
+            { status; mstate; sstate }
+        | Some z ->
+            let status = Running { r with mode = Discrete; now } in
+            let mstate = model.zset mstate z in
+            { status; mstate; sstate }
+        in
+        let fout t =
+          let t = r.now +. t in
+          model.fout mstate (r.input.u t) (f t) in
+        let out =
+          { start = r.input.start +. r.now; length = h -. r.now; u = fout } in
+        Some out, state
+  in
+  let reset _p _ = s in
+  DNode { s; step; reset }
+
diff --git a/src/lib/hsim/solver.ml b/src/lib/hsim/solver.ml
new file mode 100644
index 0000000..97b4ec6
--- /dev/null
+++ b/src/lib/hsim/solver.ml
@@ -0,0 +1,15 @@
+
+open Types
+
+(** 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 s = csolver.s, zsolver.s in
+  let step (cs, zs) h =
+    let (h', f), cs' = csolver.step cs h in
+    let (h', z), zs' = zsolver.step zs (h', f) in
+    (h', f, z), (cs', zs') in
+  let reset (ivp, zc) (cs, zs) = csolver.reset ivp cs, zsolver.reset zc zs in
+  DNode { s; step; reset }
+
diff --git a/src/lib/hsim/types.mli b/src/lib/hsim/types.mli
new file mode 100644
index 0000000..0060ef1
--- /dev/null
+++ b/src/lib/hsim/types.mli
@@ -0,0 +1,80 @@
+
+type time = float
+
+(** Input and output values are functions defined on intervals. *)
+type 'a value =
+  { start : time;
+    length : time;                     (* Relative: [end = start + length].   *)
+    u : time -> 'a }                   (* Defined on [[start, end]].          *)
+
+(** A time signal is a sequence of possibly absent α-values
+    [{ start; length; u }] where:
+  - [horizon] is a positive (possibly null) floating-point number;
+  - [u: [0, length] -> α] *)
+type 'a signal = 'a value option
+
+(** A discrete node. *)
+type ('p, 'a, 'b) dnode =
+  DNode :
+    { s : 'ds;
+      step : 'ds -> 'a -> 'b * 'ds;
+      reset : 'p -> 'ds -> 'ds;
+    } -> ('p, 'a, 'b) dnode
+
+(** A continuous node. *)
+type ('a, 'b, 'y, 'yder) cnode =
+  CNode :
+    { lsty : 'y;
+      fder : 'a -> 'y -> 'yder;
+      fout : 'a -> 'y -> 'b;
+    } -> ('a, 'b, 'y, 'yder) cnode
+
+(** A hybrid node. *)
+type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
+  HNode :
+    { s : '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.           *)
+    } -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
+
+(** An Initial Value Problem. *)
+type ('y, 'yder) ivp =
+  { init : 'y;                         (** [y₀]: initial value of y.          *)
+    fder : time -> 'y -> 'yder;        (** [dy/dt]: derivative of y.          *)
+    stop : time }                      (** Stop time.                         *)
+
+(** A zero-crossing expression. *)
+type ('y, 'zout) zc =
+  { yc : 'y;                           (** Value to watch for zero-crossings. *)
+    fzer : time -> 'y -> 'zout }       (** Zero-crossing function.            *)
+
+(** An ODE solver is a synchronous function with:
+  - 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
+
+(** 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
+
+(** 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,
+   time,
+   time * (time -> 'y) * 'zin option) dnode
diff --git a/src/lib/hsim/zls.mli b/src/lib/hsim/zls.mli
new file mode 100644
index 0000000..f28ca69
--- /dev/null
+++ b/src/lib/hsim/zls.mli
@@ -0,0 +1,7 @@
+(* This code was originally written by Timothy Bourke *)
+(* and is part of the Zelus standard library *)
+
+open Bigarray
+
+type carray = (float, float64_elt, c_layout) Array1.t
+type zarray = (int32, int32_elt, c_layout) Array1.t
diff --git a/test/dune b/test/dune
deleted file mode 100644
index 78e8da2..0000000
--- a/test/dune
+++ /dev/null
@@ -1,2 +0,0 @@
-(test
- (name test_hsim))
diff --git a/test/test_hsim.ml b/test/test_hsim.ml
deleted file mode 100644
index e69de29..0000000