diff --git a/russell_ode/src/erk_dense_out.rs b/russell_ode/src/erk_dense_out.rs
index 00f042ba..3cd6895f 100644
--- a/russell_ode/src/erk_dense_out.rs
+++ b/russell_ode/src/erk_dense_out.rs
@@ -60,7 +60,7 @@ impl ErkDenseOut {
     /// Updates the data and returns the number of function evaluations
     pub(crate) fn update<'a, F, J, A>(
         &mut self,
-        system: &mut System<'a, F, J, A>,
+        system: &System<F, J, A>,
         x: f64,
         y: &Vector,
         h: f64,
@@ -69,8 +69,8 @@ impl ErkDenseOut {
         args: &mut A,
     ) -> Result<usize, StrError>
     where
-        F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-        J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+        F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+        J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
     {
         let mut n_function_eval = 0;
 
diff --git a/russell_ode/src/euler_backward.rs b/russell_ode/src/euler_backward.rs
index b1bdabc7..b59ac09b 100644
--- a/russell_ode/src/euler_backward.rs
+++ b/russell_ode/src/euler_backward.rs
@@ -6,14 +6,14 @@ use russell_sparse::{CooMatrix, Genie, LinSolver, SparseMatrix};
 /// Implements the backward Euler (implicit) solver
 pub(crate) struct EulerBackward<'a, F, J, A>
 where
-    F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+    F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+    J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
 {
     /// Holds the parameters
     params: Params,
 
     /// ODE system
-    system: System<'a, F, J, A>,
+    system: &'a System<F, J, A>,
 
     /// Vector holding the function evaluation
     ///
@@ -38,11 +38,11 @@ where
 
 impl<'a, F, J, A> EulerBackward<'a, F, J, A>
 where
-    F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+    F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+    J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
 {
     /// Allocates a new instance
-    pub fn new(params: Params, system: System<'a, F, J, A>) -> Self {
+    pub fn new(params: Params, system: &'a System<F, J, A>) -> Self {
         let ndim = system.ndim;
         let jac_nnz = if params.newton.use_numerical_jacobian {
             ndim * ndim
@@ -67,8 +67,8 @@ where
 
 impl<'a, F, J, A> OdeSolverTrait<A> for EulerBackward<'a, F, J, A>
 where
-    F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+    F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+    J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
 {
     /// Enables dense output
     fn enable_dense_output(&mut self) -> Result<(), StrError> {
@@ -119,7 +119,9 @@ where
                 let kk = self.kk.get_coo_mut().unwrap();
                 if self.params.newton.use_numerical_jacobian || !self.system.jac_available {
                     work.bench.n_function += ndim;
-                    self.system.numerical_jacobian(kk, x_new, y_new, &self.k, h, args)?;
+                    let aux = &mut self.dy; // using dy as a workspace
+                    self.system
+                        .numerical_jacobian(kk, x_new, y_new, &self.k, h, args, aux)?;
                 } else {
                     (self.system.jacobian)(kk, x_new, y_new, h, args)?;
                 }
@@ -262,12 +264,12 @@ mod tests {
 
         // problem
         let (system, data, mut args) = Samples::kreyszig_ex4_page920();
-        let mut yfx = data.y_analytical.unwrap();
+        let yfx = data.y_analytical.unwrap();
         let ndim = system.ndim;
 
         // allocate structs
         let params = Params::new(Method::BwEuler);
-        let mut solver = EulerBackward::new(params, system);
+        let mut solver = EulerBackward::new(params, &system);
         let mut work = Workspace::new(Method::BwEuler);
 
         // check dense output availability
@@ -317,13 +319,13 @@ mod tests {
 
         // problem
         let (system, data, mut args) = Samples::kreyszig_ex4_page920();
-        let mut yfx = data.y_analytical.unwrap();
+        let yfx = data.y_analytical.unwrap();
         let ndim = system.ndim;
 
         // allocate structs
         let mut params = Params::new(Method::BwEuler);
         params.newton.use_numerical_jacobian = true;
-        let mut solver = EulerBackward::new(params, system);
+        let mut solver = EulerBackward::new(params, &system);
         let mut work = Workspace::new(Method::BwEuler);
 
         // numerical approximation
@@ -366,7 +368,7 @@ mod tests {
         let mut params = Params::new(Method::BwEuler);
         let (system, data, mut args) = Samples::kreyszig_ex4_page920();
         params.newton.n_iteration_max = 0;
-        let mut solver = EulerBackward::new(params, system);
+        let mut solver = EulerBackward::new(params, &system);
         let mut work = Workspace::new(Method::BwEuler);
         assert_eq!(
             solver.step(&mut work, data.x0, &data.y0, 0.1, &mut args).err(),
diff --git a/russell_ode/src/euler_forward.rs b/russell_ode/src/euler_forward.rs
index 4d4d34ad..a2d35522 100644
--- a/russell_ode/src/euler_forward.rs
+++ b/russell_ode/src/euler_forward.rs
@@ -5,11 +5,11 @@ use russell_sparse::CooMatrix;
 
 pub(crate) struct EulerForward<'a, F, J, A>
 where
-    F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+    F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+    J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
 {
     /// ODE system
-    system: System<'a, F, J, A>,
+    system: &'a System<F, J, A>,
 
     /// Vector holding the function evaluation
     ///
@@ -22,11 +22,11 @@ where
 
 impl<'a, F, J, A> EulerForward<'a, F, J, A>
 where
-    F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+    F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+    J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
 {
     /// Allocates a new instance
-    pub fn new(system: System<'a, F, J, A>) -> Self {
+    pub fn new(system: &'a System<F, J, A>) -> Self {
         let ndim = system.ndim;
         EulerForward {
             system,
@@ -38,8 +38,8 @@ where
 
 impl<'a, F, J, A> OdeSolverTrait<A> for EulerForward<'a, F, J, A>
 where
-    F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+    F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+    J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
 {
     /// Enables dense output
     fn enable_dense_output(&mut self) -> Result<(), StrError> {
@@ -89,11 +89,11 @@ mod tests {
 
         // problem
         let (system, data, mut args) = Samples::kreyszig_eq6_page902();
-        let mut yfx = data.y_analytical.unwrap();
+        let yfx = data.y_analytical.unwrap();
         let ndim = system.ndim;
 
         // allocate structs
-        let mut solver = EulerForward::new(system);
+        let mut solver = EulerForward::new(&system);
         let mut work = Workspace::new(Method::FwEuler);
 
         // check dense output availability
diff --git a/russell_ode/src/explicit_runge_kutta.rs b/russell_ode/src/explicit_runge_kutta.rs
index 3574c73d..69565f98 100644
--- a/russell_ode/src/explicit_runge_kutta.rs
+++ b/russell_ode/src/explicit_runge_kutta.rs
@@ -6,14 +6,14 @@ use russell_sparse::CooMatrix;
 
 pub(crate) struct ExplicitRungeKutta<'a, F, J, A>
 where
-    F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+    F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+    J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
 {
     /// Holds the parameters
     params: Params,
 
     /// ODE system
-    system: System<'a, F, J, A>,
+    system: &'a System<F, J, A>,
 
     /// Information such as implicit, embedded, etc.
     info: Information,
@@ -65,11 +65,11 @@ where
 
 impl<'a, F, J, A> ExplicitRungeKutta<'a, F, J, A>
 where
-    F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+    F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+    J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
 {
     /// Allocates a new instance
-    pub fn new(params: Params, system: System<'a, F, J, A>) -> Result<Self, StrError> {
+    pub fn new(params: Params, system: &'a System<F, J, A>) -> Result<Self, StrError> {
         // Runge-Kutta coefficients
         #[rustfmt::skip]
         let (aa, bb, cc) = match params.method {
@@ -149,8 +149,8 @@ where
 
 impl<'a, F, J, A> OdeSolverTrait<A> for ExplicitRungeKutta<'a, F, J, A>
 where
-    F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+    F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+    J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
 {
     /// Enables dense output
     fn enable_dense_output(&mut self) -> Result<(), StrError> {
@@ -385,7 +385,7 @@ mod tests {
                 None,
                 None,
             );
-            let erk = ExplicitRungeKutta::new(params, system).unwrap();
+            let erk = ExplicitRungeKutta::new(params, &system).unwrap();
             let nstage = erk.nstage;
             assert_eq!(erk.aa.dims(), (nstage, nstage));
             assert_eq!(erk.bb.dim(), nstage);
@@ -506,12 +506,12 @@ mod tests {
 
         // problem
         let (system, data, mut args) = Samples::kreyszig_eq6_page902();
-        let mut yfx = data.y_analytical.unwrap();
+        let yfx = data.y_analytical.unwrap();
         let ndim = system.ndim;
 
         // allocate structs
         let params = Params::new(Method::MdEuler); // aka the Improved Euler in Kreyszig's book
-        let mut solver = ExplicitRungeKutta::new(params, system).unwrap();
+        let mut solver = ExplicitRungeKutta::new(params, &system).unwrap();
         let mut work = Workspace::new(Method::FwEuler);
 
         // check dense output availability
@@ -589,12 +589,12 @@ mod tests {
 
         // problem
         let (system, data, mut args) = Samples::kreyszig_eq6_page902();
-        let mut yfx = data.y_analytical.unwrap();
+        let yfx = data.y_analytical.unwrap();
         let ndim = system.ndim;
 
         // allocate structs
         let params = Params::new(Method::Rk4); // aka the Classical RK in Kreyszig's book
-        let mut solver = ExplicitRungeKutta::new(params, system).unwrap();
+        let mut solver = ExplicitRungeKutta::new(params, &system).unwrap();
         let mut work = Workspace::new(Method::FwEuler);
 
         // check dense output availability
@@ -701,7 +701,7 @@ mod tests {
 
         // allocate solver
         let params = Params::new(Method::Fehlberg4);
-        let mut solver = ExplicitRungeKutta::new(params, system).unwrap();
+        let mut solver = ExplicitRungeKutta::new(params, &system).unwrap();
         let mut work = Workspace::new(Method::FwEuler);
 
         // perform one step (compute k)
diff --git a/russell_ode/src/ode_solver.rs b/russell_ode/src/ode_solver.rs
index d48a88f7..95b29ec8 100644
--- a/russell_ode/src/ode_solver.rs
+++ b/russell_ode/src/ode_solver.rs
@@ -49,10 +49,10 @@ impl<'a, A> OdeSolver<'a, A> {
     /// # Generics
     ///
     /// See [System] for an explanation of the generic parameters.
-    pub fn new<F, J>(params: Params, system: System<'a, F, J, A>) -> Result<Self, StrError>
+    pub fn new<F, J>(params: Params, system: &'a System<F, J, A>) -> Result<Self, StrError>
     where
-        F: 'a + Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-        J: 'a + Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+        F: 'a + Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+        J: 'a + Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
         A: 'a,
     {
         params.validate()?;
@@ -97,7 +97,7 @@ impl<'a, A> OdeSolver<'a, A> {
         x0: f64,
         x1: f64,
         h_equal: Option<f64>,
-        mut output: Option<&mut Output<'a>>,
+        mut output: Option<&mut Output>,
         args: &mut A,
     ) -> Result<(), StrError> {
         // check data
@@ -284,15 +284,16 @@ mod tests {
     #[test]
     fn solve_with_step_output_works() {
         // ODE system
-        let (system, data, mut args) = Samples::simple_constant();
+        let (system, data, mut args) = Samples::simple_equation_constant();
 
         // output
         let mut out = Output::new();
-        out.enable_step(&[0]).set_analytical(data.y_analytical.unwrap());
+        out.y_analytical = data.y_analytical;
+        out.enable_step(&[0]);
 
         // params and solver
         let params = Params::new(Method::FwEuler);
-        let mut solver = OdeSolver::new(params, system).unwrap();
+        let mut solver = OdeSolver::new(params, &system).unwrap();
 
         // solve the ODE system (will run with N_EQUAL_STEPS)
         let mut y = data.y0.clone();
@@ -328,7 +329,7 @@ mod tests {
     #[test]
     fn solve_with_h_equal_works() {
         // ODE system
-        let (system, mut data, mut args) = Samples::simple_constant();
+        let (system, mut data, mut args) = Samples::simple_equation_constant();
 
         // output
         let mut out = Output::new();
@@ -336,7 +337,7 @@ mod tests {
 
         // params and solver
         let params = Params::new(Method::FwEuler);
-        let mut solver = OdeSolver::new(params, system).unwrap();
+        let mut solver = OdeSolver::new(params, &system).unwrap();
         let x1 = 1.2; // => will generate 4 steps
 
         // capture error
@@ -368,16 +369,17 @@ mod tests {
     #[test]
     fn solve_with_variable_steps_works() {
         // ODE system
-        let (system, mut data, mut args) = Samples::simple_constant();
+        let (system, mut data, mut args) = Samples::simple_equation_constant();
 
         // output
         let mut out = Output::new();
-        out.enable_step(&[0]).set_analytical(data.y_analytical.unwrap());
+        out.y_analytical = data.y_analytical;
+        out.enable_step(&[0]);
 
         // params and solver
         let mut params = Params::new(Method::MdEuler);
         params.step.h_ini = 0.1;
-        let mut solver = OdeSolver::new(params, system).unwrap();
+        let mut solver = OdeSolver::new(params, &system).unwrap();
 
         // solve the ODE system
         solver
@@ -395,7 +397,7 @@ mod tests {
     #[test]
     fn solve_with_dense_output_works() {
         // ODE system
-        let (system, data, mut args) = Samples::simple_constant();
+        let (system, data, mut args) = Samples::simple_equation_constant();
 
         // output
         let mut out = Output::new();
@@ -404,7 +406,7 @@ mod tests {
         // params and solver
         let mut params = Params::new(Method::DoPri5);
         params.step.h_ini = 0.1;
-        let mut solver = OdeSolver::new(params, system).unwrap();
+        let mut solver = OdeSolver::new(params, &system).unwrap();
 
         // solve the ODE system
         let mut y = data.y0.clone();
@@ -434,7 +436,7 @@ mod tests {
     #[test]
     fn solve_handles_zero_stepsize() {
         // ODE system
-        let (system, mut data, mut args) = Samples::simple_constant();
+        let (system, mut data, mut args) = Samples::simple_equation_constant();
 
         // output
         let mut out = Output::new();
@@ -443,7 +445,7 @@ mod tests {
         // params and solver
         let mut params = Params::new(Method::DoPri5);
         params.step.h_ini = 20.0; // since the problem is linear, this stepsize will lead to a jump from x=0.0 to 1.0
-        let mut solver = OdeSolver::new(params, system).unwrap();
+        let mut solver = OdeSolver::new(params, &system).unwrap();
 
         // solve the ODE system
         solver
@@ -459,9 +461,9 @@ mod tests {
 
     #[test]
     fn solve_and_output_handle_errors() {
-        let (system, mut data, mut args) = Samples::simple_constant();
+        let (system, mut data, mut args) = Samples::simple_equation_constant();
         let params = Params::new(Method::FwEuler);
-        let mut solver = OdeSolver::new(params, system).unwrap();
+        let mut solver = OdeSolver::new(params, &system).unwrap();
         let mut out = Output::new();
         assert_eq!(out.enable_dense(-0.1, &[0]).err(), Some("h_out must be positive"));
         out.enable_dense(0.1, &[0]).unwrap();
diff --git a/russell_ode/src/output.rs b/russell_ode/src/output.rs
index 094345cd..3334bea5 100644
--- a/russell_ode/src/output.rs
+++ b/russell_ode/src/output.rs
@@ -4,7 +4,7 @@ use russell_lab::{vec_max_abs_diff, Vector};
 use std::collections::HashMap;
 
 /// Holds the (x,y) results at accepted steps or interpolated within a "dense" sequence
-pub struct Output<'a> {
+pub struct Output {
     /// Indicates whether the accepted step output is to be performed or not
     save_step: bool,
 
@@ -46,11 +46,11 @@ pub struct Output<'a> {
     /// Holds an auxiliary y vector (e.g., to compute the analytical solution or the dense output)
     y_aux: Vector,
 
-    /// Implements the analytical solution `y(x)` function
-    y_analytical: Option<Box<dyn 'a + FnMut(&mut Vector, f64)>>,
+    /// Holds a function to compute the analytical solution y(x)
+    pub y_analytical: Option<fn(&mut Vector, f64)>,
 }
 
-impl<'a> Output<'a> {
+impl Output {
     /// Allocates a new instance
     pub fn new() -> Self {
         const EMPTY: usize = 0;
@@ -130,21 +130,6 @@ impl<'a> Output<'a> {
         self.dense_h.is_some()
     }
 
-    /// Sets the analytical solution to compute the global error
-    ///
-    /// The results will be saved in the `step_global_error`
-    ///
-    /// # Input
-    ///
-    /// * `analytical` -- is a function(y, x) that computes y(x)
-    pub fn set_analytical<F>(&mut self, y_analytical: F) -> &mut Self
-    where
-        F: 'a + FnMut(&mut Vector, f64),
-    {
-        self.y_analytical = Some(Box::new(y_analytical));
-        self
-    }
-
     /// Clears all resulting arrays
     pub fn clear(&mut self) {
         self.step_h.clear();
@@ -161,7 +146,7 @@ impl<'a> Output<'a> {
     }
 
     /// Appends the results after an accepted step is computed
-    pub(crate) fn push<A>(
+    pub(crate) fn push<'a, A>(
         &mut self,
         work: &Workspace,
         x: f64,
diff --git a/russell_ode/src/radau5.rs b/russell_ode/src/radau5.rs
index 49a393e0..133c3e1a 100644
--- a/russell_ode/src/radau5.rs
+++ b/russell_ode/src/radau5.rs
@@ -9,14 +9,14 @@ use std::thread;
 /// Implements the Radau5 method
 pub(crate) struct Radau5<'a, F, J, A>
 where
-    F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+    F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+    J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
 {
     /// Holds the parameters
     params: Params,
 
     /// ODE system
-    system: System<'a, F, J, A>,
+    system: &'a System<F, J, A>,
 
     /// Holds the Jacobian matrix. J = df/dy
     jj: SparseMatrix,
@@ -104,15 +104,15 @@ where
 
 impl<'a, F, J, A> Radau5<'a, F, J, A>
 where
-    F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+    F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+    J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
 {
     /// Allocates a new instance
-    pub fn new(params: Params, system: System<'a, F, J, A>) -> Self {
+    pub fn new(params: Params, system: &'a System<F, J, A>) -> Self {
         let ndim = system.ndim;
         let symmetry = Some(system.jac_symmetry);
         let one_based = params.newton.genie == Genie::Mumps;
-        let mass_nnz = match system.mass_matrix {
+        let mass_nnz = match system.mass_matrix.as_ref() {
             Some(mass) => mass.get_info().2,
             None => ndim,
         };
@@ -186,9 +186,10 @@ where
             if self.params.newton.use_numerical_jacobian || !self.system.jac_available {
                 work.bench.n_function += self.system.ndim;
                 let y_mut = &mut self.w0; // using w[0] as a workspace
+                let aux = &mut self.dw0; // using dw0 as a workspace
                 vec_copy(y_mut, y).unwrap();
                 self.system
-                    .numerical_jacobian(jj, x, y_mut, &self.k_accepted, 1.0, args)?;
+                    .numerical_jacobian(jj, x, y_mut, &self.k_accepted, 1.0, args, aux)?;
             } else {
                 (self.system.jacobian)(jj, x, y, 1.0, args)?;
             }
@@ -202,7 +203,7 @@ where
         let gamma = GAMMA / h;
         kk_real.assign(-1.0, jj).unwrap(); // K_real = -J
         kk_comp.assign_real(-1.0, 0.0, jj).unwrap(); // K_comp = -J
-        match self.system.mass_matrix {
+        match self.system.mass_matrix.as_ref() {
             Some(mass) => {
                 kk_real.augment(gamma, mass).unwrap(); // K_real += γ M
                 kk_comp.augment_real(alpha, beta, mass).unwrap(); // K_comp += (α + βi) M
@@ -299,8 +300,8 @@ where
 
 impl<'a, F, J, A> OdeSolverTrait<A> for Radau5<'a, F, J, A>
 where
-    F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+    F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+    J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
 {
     /// Enables dense output
     fn enable_dense_output(&mut self) -> Result<(), StrError> {
@@ -394,7 +395,7 @@ where
             (self.system.function)(&mut self.k2, u2, &self.v2, args)?;
 
             // compute the right-hand side vectors
-            let (l0, l1, l2) = match self.system.mass_matrix {
+            let (l0, l1, l2) = match self.system.mass_matrix.as_ref() {
                 Some(mass) => {
                     mass.mat_vec_mul(&mut self.dw0, 1.0, &self.w0).unwrap(); // dw0 := M ⋅ w0
                     mass.mat_vec_mul(&mut self.dw1, 1.0, &self.w1).unwrap(); // dw1 := M ⋅ w1
@@ -515,7 +516,7 @@ where
         let err = &mut self.dw0; // error variable
 
         // compute ez, mez and rhs
-        match self.system.mass_matrix {
+        match self.system.mass_matrix.as_ref() {
             Some(mass) => {
                 for m in 0..ndim {
                     ez[m] = E0 * self.z0[m] + E1 * self.z1[m] + E2 * self.z2[m];
@@ -709,12 +710,12 @@ mod tests {
 
         // problem
         let (system, data, mut args) = Samples::kreyszig_ex4_page920();
-        let mut yfx = data.y_analytical.unwrap();
+        let yfx = data.y_analytical.unwrap();
         let ndim = system.ndim;
 
         // allocate structs
         let params = Params::new(Method::Radau5);
-        let mut solver = Radau5::new(params, system);
+        let mut solver = Radau5::new(params, &system);
         let mut work = Workspace::new(Method::FwEuler);
 
         // message
@@ -778,13 +779,13 @@ mod tests {
 
         // problem
         let (system, data, mut args) = Samples::kreyszig_ex4_page920();
-        let mut yfx = data.y_analytical.unwrap();
+        let yfx = data.y_analytical.unwrap();
         let ndim = system.ndim;
 
         // allocate structs
         let mut params = Params::new(Method::Radau5);
         params.newton.use_numerical_jacobian = true;
-        let mut solver = Radau5::new(params, system);
+        let mut solver = Radau5::new(params, &system);
         let mut work = Workspace::new(Method::FwEuler);
 
         // message
diff --git a/russell_ode/src/samples.rs b/russell_ode/src/samples.rs
index c9cc17a2..aa09879d 100644
--- a/russell_ode/src/samples.rs
+++ b/russell_ode/src/samples.rs
@@ -2,10 +2,10 @@ use crate::StrError;
 use crate::{HasJacobian, NoArgs, System};
 use russell_lab::math::PI;
 use russell_lab::Vector;
-use russell_sparse::CooMatrix;
+use russell_sparse::{CooMatrix, Symmetry};
 
 /// Holds data corresponding to a sample ODE problem
-pub struct SampleData<'a> {
+pub struct SampleData {
     /// Holds the initial x
     pub x0: f64,
 
@@ -18,8 +18,8 @@ pub struct SampleData<'a> {
     /// Holds the stepsize for simulations with equal-steps
     pub h_equal: Option<f64>,
 
-    /// Holds the analytical solution `y(x)`
-    pub y_analytical: Option<Box<dyn 'a + FnMut(&mut Vector, f64)>>,
+    /// Holds a function to compute the analytical solution y(x)
+    pub y_analytical: Option<fn(&mut Vector, f64)>,
 }
 
 /// Holds a collection of sample ODE problems
@@ -37,7 +37,7 @@ pub struct SampleData<'a> {
 pub struct Samples {}
 
 impl Samples {
-    /// Implements a simple ODE with a constant derivative
+    /// Implements a simple ODE with a single equation and constant derivative
     ///
     /// ```text
     /// dy
@@ -55,14 +55,13 @@ impl Samples {
     ///     * `A` -- is `NoArgs`
     /// * `data: SampleData` -- holds the initial values and the analytical solution
     /// * `args: NoArgs` -- is a placeholder variable with the arguments to F and J
-    pub fn simple_constant<'a>() -> (
+    pub fn simple_equation_constant() -> (
         System<
-            'a,
-            impl FnMut(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
-            impl FnMut(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
             NoArgs,
         >,
-        SampleData<'a>,
+        SampleData,
         NoArgs,
     ) {
         let ndim = 1;
@@ -87,9 +86,134 @@ impl Samples {
             y0: Vector::from(&[0.0]),
             x1: 1.0,
             h_equal: Some(0.2),
-            y_analytical: Some(Box::new(|y, x| {
+            y_analytical: Some(|y, x| {
                 y[0] = x;
-            })),
+            }),
+        };
+        (system, data, 0)
+    }
+
+    /// Returns a simple system with a mass matrix
+    ///
+    /// The system is:
+    ///
+    /// ```text
+    /// y0' + y1'     = -y0 + y1
+    /// y0' - y1'     =  y0 + y1
+    ///           y2' = 1/(1 + x)
+    ///
+    /// y0(0) = 1,  y1(0) = 0,  y2(0) = 0
+    /// ```
+    ///
+    /// Thus:
+    ///
+    /// ```text
+    /// M y' = f(x, y)
+    ///
+    /// with:
+    ///
+    ///     ┌          ┐       ┌           ┐
+    ///     │  1  1  0 │       │ -y0 + y1  │
+    /// M = │  1 -1  0 │   f = │  y0 + y1  │
+    ///     │  0  0  1 │       │ 1/(1 + x) │
+    ///     └          ┘       └           ┘
+    /// ```
+    ///
+    /// The Jacobian matrix is:
+    ///
+    /// ```text
+    ///          ┌          ┐
+    ///     df   │ -1  1  0 │
+    /// J = —— = │  1  1  0 │
+    ///     dy   │  0  0  0 │
+    ///          └          ┘
+    /// ```
+    ///
+    /// The analytical solution is:
+    ///
+    /// ```text
+    /// y0(x) = cos(x)
+    /// y1(x) = -sin(x)
+    /// y2(x) = log(1 + x)
+    /// ```
+    ///
+    /// # Input
+    ///
+    /// * `lower_triangle` -- Considers the symmetry of the Jacobian and Mass matrices, and generates a lower triangular representation
+    ///
+    /// # Output
+    ///
+    /// Returns `(system, data, args)` where:
+    ///
+    /// * `system: System<F, J, A>` with:
+    ///     * `F` -- is a function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)`
+    ///     * `J` -- is a function to compute the Jacobian: `(jj: &mut CooMatrix, x: f64, y: &Vector, multiplier: f64, args: &mut A)`
+    ///     * `A` -- is `NoArgs`
+    /// * `data: SampleData` -- holds the initial values
+    /// * `args: NoArgs` -- is a placeholder variable with the arguments to F and J
+    pub fn simple_system_with_mass_matrix(
+        lower_triangle: bool,
+    ) -> (
+        System<
+            impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
+            NoArgs,
+        >,
+        SampleData,
+        NoArgs,
+    ) {
+        // selected symmetry option (for both Mass and Jacobian matrices)
+        let symmetry = if lower_triangle {
+            Some(Symmetry::new_general_lower())
+        } else {
+            None
+        };
+
+        // initial values
+        let x0 = 0.0;
+        let y0 = Vector::from(&[1.0, 0.0, 0.0]);
+        let x1 = 20.0;
+
+        // ODE system
+        let ndim = 3;
+        let jac_nnz = 4;
+        let mut system = System::new(
+            ndim,
+            move |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| {
+                f[0] = -y[0] + y[1];
+                f[1] = y[0] + y[1];
+                f[2] = 1.0 / (1.0 + x);
+                Ok(())
+            },
+            move |jj: &mut CooMatrix, _x: f64, _y: &Vector, m: f64, _args: &mut NoArgs| {
+                jj.reset();
+                jj.put(0, 0, m * (-1.0)).unwrap();
+                jj.put(0, 1, m * (1.0)).unwrap();
+                jj.put(1, 0, m * (1.0)).unwrap();
+                jj.put(1, 1, m * (1.0)).unwrap();
+                Ok(())
+            },
+            HasJacobian::Yes,
+            Some(jac_nnz),
+            symmetry,
+        );
+
+        // mass matrix
+        let mass_nnz = 5;
+        system.init_mass_matrix(mass_nnz, false).unwrap();
+        system.mass_put(0, 0, 1.0).unwrap();
+        system.mass_put(0, 1, 1.0).unwrap();
+        system.mass_put(1, 0, 1.0).unwrap();
+        system.mass_put(1, 1, -1.0).unwrap();
+        system.mass_put(2, 2, 1.0).unwrap();
+
+        // control
+        let data = SampleData {
+            x0,
+            y0,
+            x1,
+            h_equal: None,
+            y_analytical: None,
         };
         (system, data, 0)
     }
@@ -116,14 +240,13 @@ impl Samples {
     ///
     /// * Kreyszig, E (2011) Advanced engineering mathematics; in collaboration with Kreyszig H,
     ///    Edward JN 10th ed 2011, Hoboken, New Jersey, Wiley
-    pub fn kreyszig_eq6_page902<'a>() -> (
+    pub fn kreyszig_eq6_page902() -> (
         System<
-            'a,
-            impl FnMut(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
-            impl FnMut(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
             NoArgs,
         >,
-        SampleData<'a>,
+        SampleData,
         NoArgs,
     ) {
         let ndim = 1;
@@ -148,9 +271,9 @@ impl Samples {
             y0: Vector::from(&[0.0]),
             x1: 1.0,
             h_equal: Some(0.2),
-            y_analytical: Some(Box::new(|y, x| {
+            y_analytical: Some(|y, x| {
                 y[0] = f64::exp(x) - x - 1.0;
-            })),
+            }),
         };
         (system, data, 0)
     }
@@ -188,14 +311,13 @@ impl Samples {
     ///
     /// * Kreyszig, E (2011) Advanced engineering mathematics; in collaboration with Kreyszig H,
     ///    Edward JN 10th ed 2011, Hoboken, New Jersey, Wiley
-    pub fn kreyszig_ex4_page920<'a>() -> (
+    pub fn kreyszig_ex4_page920() -> (
         System<
-            'a,
-            impl FnMut(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
-            impl FnMut(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
             NoArgs,
         >,
-        SampleData<'a>,
+        SampleData,
         NoArgs,
     ) {
         let ndim = 2;
@@ -223,10 +345,10 @@ impl Samples {
             y0: Vector::from(&[2.0, -10.0]),
             x1: 1.0,
             h_equal: Some(0.2),
-            y_analytical: Some(Box::new(|y, x| {
+            y_analytical: Some(|y, x| {
                 y[0] = f64::exp(-x) + f64::exp(-10.0 * x) + x;
                 y[1] = -f64::exp(-x) - 10.0 * f64::exp(-10.0 * x) + 1.0;
-            })),
+            }),
         };
         (system, data, 0)
     }
@@ -249,14 +371,13 @@ impl Samples {
     /// * Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II.
     ///   Stiff and Differential-Algebraic Problems. Second Revised Edition.
     ///   Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p
-    pub fn hairer_wanner_eq1<'a>() -> (
+    pub fn hairer_wanner_eq1() -> (
         System<
-            'a,
-            impl FnMut(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
-            impl FnMut(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
             NoArgs,
         >,
-        SampleData<'a>,
+        SampleData,
         NoArgs,
     ) {
         const L: f64 = -50.0; // lambda
@@ -282,9 +403,9 @@ impl Samples {
             y0: Vector::from(&[0.0]),
             x1: 1.5,
             h_equal: Some(1.875 / 50.0),
-            y_analytical: Some(Box::new(|y, x| {
+            y_analytical: Some(|y, x| {
                 y[0] = -L * (f64::sin(x) - L * f64::cos(x) + L * f64::exp(L * x)) / (L * L + 1.0);
-            })),
+            }),
         };
         (system, data, 0)
     }
@@ -307,14 +428,13 @@ impl Samples {
     /// * Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II.
     ///   Stiff and Differential-Algebraic Problems. Second Revised Edition.
     ///   Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p
-    pub fn robertson<'a>() -> (
+    pub fn robertson() -> (
         System<
-            'a,
-            impl FnMut(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
-            impl FnMut(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
             NoArgs,
         >,
-        SampleData<'a>,
+        SampleData,
         NoArgs,
     ) {
         let ndim = 3;
@@ -378,17 +498,16 @@ impl Samples {
     /// * Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II.
     ///   Stiff and Differential-Algebraic Problems. Second Revised Edition.
     ///   Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p
-    pub fn van_der_pol<'a>(
+    pub fn van_der_pol(
         epsilon: Option<f64>,
         stationary: bool,
     ) -> (
         System<
-            'a,
-            impl FnMut(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
-            impl FnMut(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
             NoArgs,
         >,
-        SampleData<'a>,
+        SampleData,
         NoArgs,
     ) {
         let mut eps = match epsilon {
@@ -477,14 +596,13 @@ impl Samples {
     /// * Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II.
     ///   Stiff and Differential-Algebraic Problems. Second Revised Edition.
     ///   Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p
-    pub fn arenstorf<'a>() -> (
+    pub fn arenstorf() -> (
         System<
-            'a,
-            impl FnMut(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
-            impl FnMut(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
             NoArgs,
         >,
-        SampleData<'a>,
+        SampleData,
         NoArgs,
     ) {
         const MU: f64 = 0.012277471;
@@ -555,7 +673,7 @@ impl Samples {
     ///
     /// # Output
     ///
-    /// Returns `(system, data, args, gen_mass_matrix)` where:
+    /// Returns `(system, data, args)` where:
     ///
     /// * `system: System<F, J, A>` with:
     ///     * `F` -- is a function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)`
@@ -563,25 +681,20 @@ impl Samples {
     ///     * `A` -- is `NoArgs`
     /// * `data: SampleData` -- holds the initial values
     /// * `args: NoArgs` -- is a placeholder variable with the arguments to F and J
-    /// * `gen_mass_matrix: fn(one_based: bool) -> CooMatrix` -- is a function to generate the mass matrix.
-    ///    Note: the mass matrix needs to be allocated externally because a reference to it is required by System.
     ///
     /// # Reference
     ///
     /// * Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II.
     ///   Stiff and Differential-Algebraic Problems. Second Revised Edition.
     ///   Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p
-
-    pub fn amplifier1t<'a>() -> (
+    pub fn amplifier1t() -> (
         System<
-            'a,
-            impl FnMut(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
-            impl FnMut(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
+            impl Fn(&mut CooMatrix, f64, &Vector, f64, &mut NoArgs) -> Result<(), StrError>,
             NoArgs,
         >,
-        SampleData<'a>,
+        SampleData,
         NoArgs,
-        fn(bool) -> CooMatrix,
     ) {
         // constants
         const ALPHA: f64 = 0.99;
@@ -602,7 +715,7 @@ impl Samples {
         // ODE system
         let ndim = 5;
         let jac_nnz = 9;
-        let system = System::new(
+        let mut system = System::new(
             ndim,
             move |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| {
                 let ue = C * f64::sin(D * x);
@@ -634,26 +747,22 @@ impl Samples {
         );
 
         // function that generates the mass matrix
-        let gen_mass_matrix = |one_based: bool| -> CooMatrix {
-            const C1: f64 = 1e-6;
-            const C2: f64 = 2e-6;
-            const C3: f64 = 3e-6;
-            let ndim = 5;
-            let nnz = 9;
-            let mut mass = CooMatrix::new(ndim, ndim, nnz, None, one_based).unwrap();
-            mass.put(0, 0, -C1).unwrap();
-            mass.put(0, 1, C1).unwrap();
-            mass.put(1, 0, C1).unwrap();
-            mass.put(1, 1, -C1).unwrap();
-            mass.put(2, 2, -C2).unwrap();
-            mass.put(3, 3, -C3).unwrap();
-            mass.put(3, 4, C3).unwrap();
-            mass.put(4, 3, C3).unwrap();
-            mass.put(4, 4, -C3).unwrap();
-            mass
-        };
+        const C1: f64 = 1e-6;
+        const C2: f64 = 2e-6;
+        const C3: f64 = 3e-6;
+        let mass_nnz = 9;
+        system.init_mass_matrix(mass_nnz, false).unwrap();
+        system.mass_put(0, 0, -C1).unwrap();
+        system.mass_put(0, 1, C1).unwrap();
+        system.mass_put(1, 0, C1).unwrap();
+        system.mass_put(1, 1, -C1).unwrap();
+        system.mass_put(2, 2, -C2).unwrap();
+        system.mass_put(3, 3, -C3).unwrap();
+        system.mass_put(3, 4, C3).unwrap();
+        system.mass_put(4, 3, C3).unwrap();
+        system.mass_put(4, 4, -C3).unwrap();
 
-        // results
+        // control
         let data = SampleData {
             x0,
             y0,
@@ -661,7 +770,7 @@ impl Samples {
             h_equal: None,
             y_analytical: None,
         };
-        (system, data, 0, gen_mass_matrix)
+        (system, data, 0)
     }
 }
 
@@ -674,9 +783,9 @@ mod tests {
     use russell_lab::{deriv_central5, mat_approx_eq, vec_approx_eq, Matrix, Vector};
     use russell_sparse::{CooMatrix, Symmetry};
 
-    fn numerical_jacobian<F>(ndim: usize, x0: f64, y0: Vector, mut function: F, multiplier: f64) -> Matrix
+    fn numerical_jacobian<F>(ndim: usize, x0: f64, y0: Vector, function: F, multiplier: f64) -> Matrix
     where
-        F: FnMut(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
+        F: Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>,
     {
         struct Extra {
             x: f64,
@@ -715,7 +824,7 @@ mod tests {
     #[test]
     fn simple_constant_works() {
         let multiplier = 2.0;
-        let (mut system, mut data, mut args) = Samples::simple_constant();
+        let (system, mut data, mut args) = Samples::simple_equation_constant();
 
         // check initial values
         if let Some(y_ana) = data.y_analytical.as_mut() {
@@ -743,7 +852,7 @@ mod tests {
     #[test]
     fn kreyszig_eq6_page902_works() {
         let multiplier = 2.0;
-        let (mut system, mut data, mut args) = Samples::kreyszig_eq6_page902();
+        let (system, mut data, mut args) = Samples::kreyszig_eq6_page902();
 
         // check initial values
         if let Some(y_ana) = data.y_analytical.as_mut() {
@@ -771,7 +880,7 @@ mod tests {
     #[test]
     fn kreyszig_ex4_page920() {
         let multiplier = 2.0;
-        let (mut system, mut data, mut args) = Samples::kreyszig_ex4_page920();
+        let (system, mut data, mut args) = Samples::kreyszig_ex4_page920();
 
         // check initial values
         if let Some(y_ana) = data.y_analytical.as_mut() {
@@ -799,7 +908,7 @@ mod tests {
     #[test]
     fn hairer_wanner_eq1_works() {
         let multiplier = 2.0;
-        let (mut system, mut data, mut args) = Samples::hairer_wanner_eq1();
+        let (system, mut data, mut args) = Samples::hairer_wanner_eq1();
 
         // check initial values
         if let Some(y_ana) = data.y_analytical.as_mut() {
@@ -827,7 +936,7 @@ mod tests {
     #[test]
     fn robertson_works() {
         let multiplier = 2.0;
-        let (mut system, data, mut args) = Samples::robertson();
+        let (system, data, mut args) = Samples::robertson();
 
         // compute the analytical Jacobian matrix
         let symmetry = Some(system.jac_symmetry);
@@ -847,7 +956,7 @@ mod tests {
     #[test]
     fn van_der_pol_works() {
         let multiplier = 2.0;
-        let (mut system, data, mut args) = Samples::van_der_pol(None, false);
+        let (system, data, mut args) = Samples::van_der_pol(None, false);
 
         // compute the analytical Jacobian matrix
         let symmetry = Some(system.jac_symmetry);
@@ -867,7 +976,7 @@ mod tests {
     #[test]
     fn van_der_pol_works_stationary() {
         let multiplier = 3.0;
-        let (mut system, data, mut args) = Samples::van_der_pol(None, true);
+        let (system, data, mut args) = Samples::van_der_pol(None, true);
 
         // compute the analytical Jacobian matrix
         let symmetry = Some(system.jac_symmetry);
@@ -887,7 +996,7 @@ mod tests {
     #[test]
     fn arenstorf_works() {
         let multiplier = 1.5;
-        let (mut system, data, mut args) = Samples::arenstorf();
+        let (system, data, mut args) = Samples::arenstorf();
 
         // compute the analytical Jacobian matrix
         let symmetry = Some(system.jac_symmetry);
@@ -907,7 +1016,7 @@ mod tests {
     #[test]
     fn amplifier1t_works() {
         let multiplier = 2.0;
-        let (mut system, data, mut args, gen_mass_matrix) = Samples::amplifier1t();
+        let (system, data, mut args) = Samples::amplifier1t();
 
         // compute the analytical Jacobian matrix
         let symmetry = Some(system.jac_symmetry);
@@ -923,8 +1032,8 @@ mod tests {
         println!("{:.15}", num);
         mat_approx_eq(&ana, &num, 1e-13);
 
-        // generate the mass matrix
-        let mass = gen_mass_matrix(false);
+        // check the mass matrix
+        let mass = system.mass_matrix.unwrap();
         println!("{}", mass.as_dense());
         let ndim = system.ndim;
         let nnz_mass = 5 + 4;
diff --git a/russell_ode/src/system.rs b/russell_ode/src/system.rs
index 66cca33f..54dcb67d 100644
--- a/russell_ode/src/system.rs
+++ b/russell_ode/src/system.rs
@@ -6,35 +6,43 @@ use std::marker::PhantomData;
 /// Indicates that the system functions do not require extra arguments
 pub type NoArgs = u8;
 
-/// Defines the system of ODEs
+/// Defines a system of first order ordinary differential equations (ODE) or a differential-algebraic system (DAE) of Index-1
 ///
 /// The system is defined by:
 ///
 /// ```text
-/// d{y}
-/// ———— = {f}(x, {y})
-///  dx
-/// where x is a scalar and {y} and {f} are vectors
+///     d{y}
+/// [M] ———— = {f}(x, {y})
+///      dx
 /// ```
 ///
+/// where `x` is the independent scalar variable (e.g., time), `{y}` is the solution vector,
+/// `{f}` is the right-hand side vector, and `[M]` is the so-called "mass matrix".
+///
+/// **Note:** The mass matrix is optional and need not be specified.
+///
 /// The Jacobian is defined by:
 ///
 /// ```text
 ///               ∂{f}
 /// [J](x, {y}) = ————
 ///               ∂{y}
-/// where [J] is the Jacobian matrix
 /// ```
 ///
+/// where `[J]` is the Jacobian matrix.
+///
 /// # Generics
 ///
-/// * `F` -- is a function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)`
-/// * `J` -- is a function to compute the Jacobian: `(jj: &mut CooMatrix, x: f64, y: &Vector, multiplier: f64, args: &mut A)`
-/// * `A` -- is a generic argument to assist in the `F` and `J` functions
+/// The generic arguments here are:
+///
+/// * `F` -- function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)`
+/// * `J` -- function to compute the Jacobian: `(jj: &mut CooMatrix, x: f64, y: &Vector, multiplier: f64, args: &mut A)`
+/// * `A` -- generic argument to assist in the `F` and `J` functions. It may be simply the [NoArgs] type indicating that no arguments are needed.
 ///
 /// # Important
 ///
-/// The `multiplier` parameter in the Jacobian function `J` must be used to scale the Jacobian matrix; e.g.,
+/// The internal implementation requires that the `multiplier` parameter in
+/// the Jacobian function `J` be used used to scale the Jacobian matrix. For example:
 ///
 /// ```text
 /// |jj: &mut CooMatrix, x: f64, y: &Vector, multiplier: f64, args: &mut Args| {
@@ -43,10 +51,10 @@ pub type NoArgs = u8;
 ///     Ok(())
 /// },
 /// ```
-pub struct System<'a, F, J, A>
+pub struct System<F, J, A>
 where
-    F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+    F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+    J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
 {
     /// System dimension
     pub(crate) ndim: usize,
@@ -67,19 +75,16 @@ where
     pub(crate) jac_symmetry: Symmetry,
 
     /// Holds the mass matrix
-    pub(crate) mass_matrix: Option<&'a CooMatrix>,
-
-    /// workspace for numerical Jacobian
-    work: Vector,
+    pub(crate) mass_matrix: Option<CooMatrix>,
 
     /// Handle generic argument
     phantom: PhantomData<fn() -> A>,
 }
 
-impl<'a, F, J, A> System<'a, F, J, A>
+impl<'a, F, J, A> System<F, J, A>
 where
-    F: Send + FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    J: Send + FnMut(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
+    F: Send + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
+    J: Send + Fn(&mut CooMatrix, f64, &Vector, f64, &mut A) -> Result<(), StrError>,
 {
     /// Allocates a new instance
     ///
@@ -102,10 +107,7 @@ where
         has_ana_jacobian: HasJacobian,
         jac_nnz: Option<usize>,
         jac_symmetry: Option<Symmetry>,
-    ) -> Self
-    where
-        F: FnMut(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    {
+    ) -> Self {
         let jac_available = match has_ana_jacobian {
             HasJacobian::Yes => true,
             HasJacobian::No => false,
@@ -118,30 +120,45 @@ where
             jac_nnz: if let Some(n) = jac_nnz { n } else { ndim * ndim },
             jac_symmetry: if let Some(s) = jac_symmetry { s } else { Symmetry::No },
             mass_matrix: None,
-            work: Vector::new(ndim),
             phantom: PhantomData,
         }
     }
 
-    /// Specifies the mass matrix
-    pub fn set_mass_matrix(&mut self, mass: &'a CooMatrix) -> Result<(), StrError> {
-        let (nrow, ncol, nnz, symmetry) = mass.get_info();
-        if nrow != ncol {
-            return Err("the mass matrix must be square");
-        }
-        if nnz < nrow {
-            return Err("the number of non-zero values in the mass matrix must be ≥ ndim");
-        }
-        if nrow != self.ndim {
-            return Err("the dimension of the mass matrix must be = ndim");
-        }
-        if symmetry != self.jac_symmetry {
-            return Err("the symmetry of the mass matrix must be = Jacobian symmetry");
-        }
-        self.mass_matrix = Some(mass);
+    /// Initializes and enables the mass matrix
+    ///
+    /// **Note:** Later, call [System::mass_put] to "put" elements in the mass matrix.
+    ///
+    /// # Input
+    ///
+    /// * `max_nnz` -- Max number of non-zero values
+    /// * `one_based` -- Flag indicating whether the Sparse matrix can be used with Fortran code (e.g., MUMPS) or not.
+    ///   Make sure that this flag is the same used for the Jacobian matrix.
+    pub fn init_mass_matrix(&mut self, max_nnz: usize, one_based: bool) -> Result<(), StrError> {
+        let sym = if self.jac_symmetry == Symmetry::No {
+            None
+        } else {
+            Some(self.jac_symmetry)
+        };
+        self.mass_matrix = Some(CooMatrix::new(self.ndim, self.ndim, max_nnz, sym, one_based).unwrap());
         Ok(())
     }
 
+    /// Puts a new element in the mass matrix (duplicates allowed)
+    ///
+    /// See also [russell_sparse::CooMatrix::put].
+    ///
+    /// # Input
+    ///
+    /// * `i` -- row index (indices start at zero; zero-based)
+    /// * `j` -- column index (indices start at zero; zero-based)
+    /// * `value` -- the value M(i,j)
+    pub fn mass_put(&mut self, i: usize, j: usize, value: f64) -> Result<(), StrError> {
+        match self.mass_matrix.as_mut() {
+            Some(mass) => mass.put(i, j, value),
+            None => Err("mass matrix has not been initialized/enabled"),
+        }
+    }
+
     /// Returns the dimension of the ODE system
     pub fn get_ndim(&self) -> usize {
         self.ndim
@@ -157,23 +174,25 @@ where
     ///
     /// **Note:** Will call `function` ndim times.
     pub(crate) fn numerical_jacobian(
-        &mut self,
+        &self,
         jj: &mut CooMatrix,
         x: f64,
         y: &mut Vector,
         fxy: &Vector,
         multiplier: f64,
         args: &mut A,
+        aux: &mut Vector,
     ) -> Result<(), StrError> {
+        assert_eq!(aux.dim(), self.ndim);
         const THRESHOLD: f64 = 1e-5;
         jj.reset();
         for j in 0..self.ndim {
             let yj_original = y[j]; // create copy
             let delta_yj = f64::sqrt(f64::EPSILON * f64::max(THRESHOLD, f64::abs(y[j])));
             y[j] += delta_yj; // y[j] := y[j] + Δy
-            (self.function)(&mut self.work, x, y, args)?; // work := f(x, y + Δy)
+            (self.function)(aux, x, y, args)?; // work := f(x, y + Δy)
             for i in 0..self.ndim {
-                let delta_fi = self.work[i] - fxy[i]; // compute Δf[..]
+                let delta_fi = aux[i] - fxy[i]; // compute Δf[..]
                 jj.put(i, j, multiplier * delta_fi / delta_yj).unwrap(); // Δfi/Δfj
             }
             y[j] = yj_original; // restore value
@@ -206,17 +225,18 @@ mod tests {
 
     #[test]
     fn ode_system_most_none_works() {
-        let mut n_function_eval = 0;
         struct Args {
+            n_function_eval: usize,
             more_data_goes_here: bool,
         }
         let mut args = Args {
+            n_function_eval: 0,
             more_data_goes_here: false,
         };
-        let mut ode = System::new(
+        let ode = System::new(
             2,
             |f, x, y, args: &mut Args| {
-                n_function_eval += 1;
+                args.n_function_eval += 1;
                 f[0] = -x * y[1];
                 f[1] = x * y[0];
                 args.more_data_goes_here = true;
@@ -240,37 +260,36 @@ mod tests {
             Err("analytical Jacobian is not available")
         );
         // check
-        println!("n_function_eval = {}", n_function_eval);
-        assert_eq!(n_function_eval, 1);
+        println!("n_function_eval = {}", args.n_function_eval);
+        assert_eq!(args.n_function_eval, 1);
         assert_eq!(args.more_data_goes_here, true);
     }
 
     #[test]
     fn ode_system_some_none_works() {
-        let mut mass = CooMatrix::new(2, 2, 2, None, false).unwrap();
-        mass.put(0, 0, 1.0).unwrap();
-        mass.put(1, 1, 1.0).unwrap();
-        let mut n_function_eval = 0;
-        let mut n_jacobian_eval = 0;
         struct Args {
+            n_function_eval: usize,
+            n_jacobian_eval: usize,
             more_data_goes_here_fn: bool,
             more_data_goes_here_jj: bool,
         }
         let mut args = Args {
+            n_function_eval: 0,
+            n_jacobian_eval: 0,
             more_data_goes_here_fn: false,
             more_data_goes_here_jj: false,
         };
         let mut ode = System::new(
             2,
             |f, x, y, args: &mut Args| {
-                n_function_eval += 1;
+                args.n_function_eval += 1;
                 f[0] = -x * y[1];
                 f[1] = x * y[0];
                 args.more_data_goes_here_fn = true;
                 Ok(())
             },
             |jj, x, _y, _multiplier, args: &mut Args| {
-                n_jacobian_eval += 1;
+                args.n_jacobian_eval += 1;
                 jj.reset();
                 jj.put(0, 1, -x).unwrap();
                 jj.put(1, 0, x).unwrap();
@@ -281,11 +300,12 @@ mod tests {
             Some(2),
             None,
         );
-        // ode.set_analytical_solution(Box::new(|y, x| {
-        //     y[0] = f64::cos(x * x / 2.0) - 2.0 * f64::sin(x * x / 2.0);
-        //     y[1] = 2.0 * f64::cos(x * x / 2.0) + f64::sin(x * x / 2.0);
-        // }));
-        ode.set_mass_matrix(&mass).unwrap();
+        // analytical_solution:
+        // y[0] = f64::cos(x * x / 2.0) - 2.0 * f64::sin(x * x / 2.0);
+        // y[1] = 2.0 * f64::cos(x * x / 2.0) + f64::sin(x * x / 2.0);
+        ode.init_mass_matrix(2, false).unwrap(); // diagonal mass matrix => OK, but not needed
+        ode.mass_put(0, 0, 1.0).unwrap();
+        ode.mass_put(1, 1, 1.0).unwrap();
         // call system function
         let x = 0.0;
         let y = Vector::new(2);
@@ -296,10 +316,10 @@ mod tests {
         let m = 1.0;
         (ode.jacobian)(&mut jj, x, &y, m, &mut args).unwrap();
         // check
-        println!("n_function_eval = {}", n_function_eval);
-        println!("n_jacobian_eval = {}", n_jacobian_eval);
-        assert_eq!(n_function_eval, 1);
-        assert_eq!(n_jacobian_eval, 1);
+        println!("n_function_eval = {}", args.n_function_eval);
+        println!("n_jacobian_eval = {}", args.n_jacobian_eval);
+        assert_eq!(args.n_function_eval, 1);
+        assert_eq!(args.n_jacobian_eval, 1);
         assert_eq!(args.more_data_goes_here_fn, true);
         assert_eq!(args.more_data_goes_here_jj, true);
     }
diff --git a/russell_ode/tests/test_bweuler.rs b/russell_ode/tests/test_bweuler.rs
index 0e549fc7..f29e2d4b 100644
--- a/russell_ode/tests/test_bweuler.rs
+++ b/russell_ode/tests/test_bweuler.rs
@@ -9,7 +9,7 @@ fn test_bweuler_hairer_wanner_eq1() {
 
     // set configuration parameters
     let params = Params::new(Method::BwEuler);
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
 
     // solve the ODE system
     solver
@@ -24,7 +24,7 @@ fn test_bweuler_hairer_wanner_eq1() {
     assert_eq!(stat.h_accepted, data.h_equal.unwrap());
 
     // compare with the analytical solution
-    let mut analytical = data.y_analytical.unwrap();
+    let analytical = data.y_analytical.unwrap();
     let mut y1_correct = Vector::new(ndim);
     analytical(&mut y1_correct, data.x1);
     approx_eq(data.y0[0], y1_correct[0], 5e-5);
@@ -53,7 +53,7 @@ fn test_bweuler_hairer_wanner_eq1_num_jac() {
     params.newton.use_numerical_jacobian = true;
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver
         .solve(&mut data.y0, data.x0, data.x1, data.h_equal, None, &mut args)
         .unwrap();
@@ -66,7 +66,7 @@ fn test_bweuler_hairer_wanner_eq1_num_jac() {
     assert_eq!(stat.h_accepted, data.h_equal.unwrap());
 
     // compare with the analytical solution
-    let mut analytical = data.y_analytical.unwrap();
+    let analytical = data.y_analytical.unwrap();
     let mut y1_correct = Vector::new(ndim);
     analytical(&mut y1_correct, data.x1);
     approx_eq(data.y0[0], y1_correct[0], 5e-5);
@@ -95,7 +95,7 @@ fn test_bweuler_hairer_wanner_eq1_modified_newton() {
     params.bweuler.use_modified_newton = true;
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver
         .solve(&mut data.y0, data.x0, data.x1, data.h_equal, None, &mut args)
         .unwrap();
@@ -108,7 +108,7 @@ fn test_bweuler_hairer_wanner_eq1_modified_newton() {
     assert_eq!(stat.h_accepted, data.h_equal.unwrap());
 
     // compare with the analytical solution
-    let mut analytical = data.y_analytical.unwrap();
+    let analytical = data.y_analytical.unwrap();
     let mut y1_correct = Vector::new(ndim);
     analytical(&mut y1_correct, data.x1);
     approx_eq(data.y0[0], y1_correct[0], 5e-5);
diff --git a/russell_ode/tests/test_dopri5_arenstorf.rs b/russell_ode/tests/test_dopri5_arenstorf.rs
index ddd677f7..b7ec14bc 100644
--- a/russell_ode/tests/test_dopri5_arenstorf.rs
+++ b/russell_ode/tests/test_dopri5_arenstorf.rs
@@ -16,7 +16,7 @@ fn test_dopri5_arenstorf() {
     out.enable_dense(1.0, &[0, 1, 2, 3]).unwrap();
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver
         .solve(&mut data.y0, data.x0, data.x1, None, Some(&mut out), &mut args)
         .unwrap();
diff --git a/russell_ode/tests/test_dopri5_arenstorf_debug.rs b/russell_ode/tests/test_dopri5_arenstorf_debug.rs
index d4d1bdc9..946c7794 100644
--- a/russell_ode/tests/test_dopri5_arenstorf_debug.rs
+++ b/russell_ode/tests/test_dopri5_arenstorf_debug.rs
@@ -13,7 +13,7 @@ fn test_dopri5_arenstorf_debug() {
     params.debug = true;
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver
         .solve(&mut data.y0, data.x0, data.x1, None, None, &mut args)
         .unwrap();
diff --git a/russell_ode/tests/test_dopri5_hairer_wanner_eq1.rs b/russell_ode/tests/test_dopri5_hairer_wanner_eq1.rs
index 19e2c56c..d9cb5433 100644
--- a/russell_ode/tests/test_dopri5_hairer_wanner_eq1.rs
+++ b/russell_ode/tests/test_dopri5_hairer_wanner_eq1.rs
@@ -16,7 +16,7 @@ fn test_dopri5_hairer_wanner_eq1() {
     out.enable_dense(0.1, &[0]).unwrap();
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver
         .solve(&mut data.y0, data.x0, data.x1, None, Some(&mut out), &mut args)
         .unwrap();
@@ -28,7 +28,7 @@ fn test_dopri5_hairer_wanner_eq1() {
     approx_eq(data.y0[0], 9.063921649310544E-02, 1e-13);
 
     // compare with the analytical solution
-    let mut analytical = data.y_analytical.unwrap();
+    let analytical = data.y_analytical.unwrap();
     let mut y1_correct = Vector::new(ndim);
     analytical(&mut y1_correct, data.x1);
     approx_eq(data.y0[0], y1_correct[0], 4e-5);
diff --git a/russell_ode/tests/test_dopri5_van_der_pol_debug.rs b/russell_ode/tests/test_dopri5_van_der_pol_debug.rs
index 1479a130..4ff4b59e 100644
--- a/russell_ode/tests/test_dopri5_van_der_pol_debug.rs
+++ b/russell_ode/tests/test_dopri5_van_der_pol_debug.rs
@@ -24,7 +24,7 @@ fn test_dopri5_van_der_pol_debug() {
     let mut y0 = Vector::from(&[2.0, 0.0]);
     let x0 = 0.0;
     let x1 = 2.0;
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver.solve(&mut y0, x0, x1, None, Some(&mut out), &mut args).unwrap();
 
     // get statistics
diff --git a/russell_ode/tests/test_dopri8_van_der_pol.rs b/russell_ode/tests/test_dopri8_van_der_pol.rs
index 6b5dc0fc..20adcb59 100644
--- a/russell_ode/tests/test_dopri8_van_der_pol.rs
+++ b/russell_ode/tests/test_dopri8_van_der_pol.rs
@@ -20,7 +20,7 @@ fn test_dopri8_van_der_pol() {
     let mut y0 = Vector::from(&[2.0, 0.0]);
     let x0 = 0.0;
     let x1 = 2.0;
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver.solve(&mut y0, x0, x1, None, Some(&mut out), &mut args).unwrap();
 
     // get statistics
diff --git a/russell_ode/tests/test_dopri8_van_der_pol_debug.rs b/russell_ode/tests/test_dopri8_van_der_pol_debug.rs
index 7339f81b..51f056fa 100644
--- a/russell_ode/tests/test_dopri8_van_der_pol_debug.rs
+++ b/russell_ode/tests/test_dopri8_van_der_pol_debug.rs
@@ -24,7 +24,7 @@ fn test_dopri8_van_der_pol_debug() {
     let mut y0 = Vector::from(&[2.0, 0.0]);
     let x0 = 0.0;
     let x1 = 2.0;
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver.solve(&mut y0, x0, x1, None, Some(&mut out), &mut args).unwrap();
 
     // get statistics
diff --git a/russell_ode/tests/test_fweuler.rs b/russell_ode/tests/test_fweuler.rs
index 8bfd8d38..41100283 100644
--- a/russell_ode/tests/test_fweuler.rs
+++ b/russell_ode/tests/test_fweuler.rs
@@ -11,7 +11,7 @@ fn test_fweuler_hairer_wanner_eq1() {
     let params = Params::new(Method::FwEuler);
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver
         .solve(&mut data.y0, data.x0, data.x1, data.h_equal, None, &mut args)
         .unwrap();
@@ -24,7 +24,7 @@ fn test_fweuler_hairer_wanner_eq1() {
     assert_eq!(stat.h_accepted, data.h_equal.unwrap());
 
     // compare with the analytical solution
-    let mut analytical = data.y_analytical.unwrap();
+    let analytical = data.y_analytical.unwrap();
     let mut y1_correct = Vector::new(ndim);
     analytical(&mut y1_correct, data.x1);
     approx_eq(data.y0[0], y1_correct[0], 0.004753);
diff --git a/russell_ode/tests/test_mdeuler.rs b/russell_ode/tests/test_mdeuler.rs
index ba7ff184..4a5b0f71 100644
--- a/russell_ode/tests/test_mdeuler.rs
+++ b/russell_ode/tests/test_mdeuler.rs
@@ -12,7 +12,7 @@ fn test_mdeuler_hairer_wanner_eq1() {
     params.step.h_ini = 1e-4;
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver
         .solve(&mut data.y0, data.x0, data.x1, None, None, &mut args)
         .unwrap();
@@ -24,7 +24,7 @@ fn test_mdeuler_hairer_wanner_eq1() {
     approx_eq(data.y0[0], 0.09062475637905158, 1e-16);
 
     // compare with the analytical solution
-    let mut analytical = data.y_analytical.unwrap();
+    let analytical = data.y_analytical.unwrap();
     let mut y1_correct = Vector::new(ndim);
     analytical(&mut y1_correct, data.x1);
     approx_eq(data.y0[0], y1_correct[0], 1e-4);
diff --git a/russell_ode/tests/test_radau5_amplifier1t.rs b/russell_ode/tests/test_radau5_amplifier1t.rs
index 560eca25..6f062639 100644
--- a/russell_ode/tests/test_radau5_amplifier1t.rs
+++ b/russell_ode/tests/test_radau5_amplifier1t.rs
@@ -4,24 +4,19 @@ use russell_ode::{Method, OdeSolver, Output, Params, Samples};
 #[test]
 fn test_radau5_amplifier1t() {
     // get get ODE system
-    let (mut system, mut data, mut args, gen_mass) = Samples::amplifier1t();
+    let (system, mut data, mut args) = Samples::amplifier1t();
 
     // set configuration parameters
     let mut params = Params::new(Method::Radau5);
     params.step.h_ini = 1e-6;
     params.set_tolerances(1e-4, 1e-4, None).unwrap();
 
-    // mass matrix
-    let one_based = false; // change to true if using MUMPS
-    let mass = gen_mass(one_based);
-    system.set_mass_matrix(&mass).unwrap();
-
     // enable output of accepted steps
     let mut out = Output::new();
     out.enable_dense(0.001, &[0, 4]).unwrap();
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver
         .solve(&mut data.y0, data.x0, data.x1, None, Some(&mut out), &mut args)
         .unwrap();
diff --git a/russell_ode/tests/test_radau5_hairer_wanner_eq1.rs b/russell_ode/tests/test_radau5_hairer_wanner_eq1.rs
index 0af24442..cd0eb3d3 100644
--- a/russell_ode/tests/test_radau5_hairer_wanner_eq1.rs
+++ b/russell_ode/tests/test_radau5_hairer_wanner_eq1.rs
@@ -16,7 +16,7 @@ fn test_radau5_hairer_wanner_eq1() {
     out.enable_dense(0.2, &[0]).unwrap();
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver
         .solve(&mut data.y0, data.x0, data.x1, None, Some(&mut out), &mut args)
         .unwrap();
@@ -29,7 +29,7 @@ fn test_radau5_hairer_wanner_eq1() {
     approx_eq(stat.h_accepted, 1.272673814374611E+00, 1e-12);
 
     // compare with the analytical solution
-    let mut analytical = data.y_analytical.unwrap();
+    let analytical = data.y_analytical.unwrap();
     let mut y1_correct = Vector::new(ndim);
     analytical(&mut y1_correct, data.x1);
     approx_eq(data.y0[0], y1_correct[0], 3e-5);
diff --git a/russell_ode/tests/test_radau5_hairer_wanner_eq1_debug.rs b/russell_ode/tests/test_radau5_hairer_wanner_eq1_debug.rs
index 3557d82a..1c5f4047 100644
--- a/russell_ode/tests/test_radau5_hairer_wanner_eq1_debug.rs
+++ b/russell_ode/tests/test_radau5_hairer_wanner_eq1_debug.rs
@@ -13,7 +13,7 @@ fn test_radau5_hairer_wanner_eq1_debug() {
     params.debug = true;
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver
         .solve(&mut data.y0, data.x0, data.x1, None, None, &mut args)
         .unwrap();
@@ -26,7 +26,7 @@ fn test_radau5_hairer_wanner_eq1_debug() {
     approx_eq(stat.h_accepted, 1.272673814374611E+00, 1e-12);
 
     // compare with the analytical solution
-    let mut analytical = data.y_analytical.unwrap();
+    let analytical = data.y_analytical.unwrap();
     let mut y1_correct = Vector::new(ndim);
     analytical(&mut y1_correct, data.x1);
     approx_eq(data.y0[0], y1_correct[0], 3e-5);
diff --git a/russell_ode/tests/test_radau5_robertson.rs b/russell_ode/tests/test_radau5_robertson.rs
index 8ff7079d..c92f98af 100644
--- a/russell_ode/tests/test_radau5_robertson.rs
+++ b/russell_ode/tests/test_radau5_robertson.rs
@@ -16,7 +16,7 @@ fn test_radau5_robertson() {
     out.enable_step(&[0, 1, 2]);
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver
         .solve(&mut data.y0, data.x0, data.x1, None, Some(&mut out), &mut args)
         .unwrap();
diff --git a/russell_ode/tests/test_radau5_robertson_debug.rs b/russell_ode/tests/test_radau5_robertson_debug.rs
index f2b4e090..8fd38350 100644
--- a/russell_ode/tests/test_radau5_robertson_debug.rs
+++ b/russell_ode/tests/test_radau5_robertson_debug.rs
@@ -13,7 +13,7 @@ fn test_radau5_robertson_debug() {
     params.debug = true;
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver
         .solve(&mut data.y0, data.x0, data.x1, None, None, &mut args)
         .unwrap();
diff --git a/russell_ode/tests/test_radau5_robertson_small_h.rs b/russell_ode/tests/test_radau5_robertson_small_h.rs
index 7282142d..59b6e972 100644
--- a/russell_ode/tests/test_radau5_robertson_small_h.rs
+++ b/russell_ode/tests/test_radau5_robertson_small_h.rs
@@ -15,7 +15,7 @@ fn test_radau5_robertson_small_h() {
     params.set_tolerances(1e-2, 1e-2, None).unwrap();
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     let res = solver.solve(&mut data.y0, data.x0, data.x1, None, None, &mut args);
     assert_eq!(res.err(), Some("the stepsize becomes too small"));
     println!("ERROR: THE STEPSIZE BECOMES TOO SMALL");
diff --git a/russell_ode/tests/test_radau5_van_der_pol.rs b/russell_ode/tests/test_radau5_van_der_pol.rs
index 43416b2a..cd70597c 100644
--- a/russell_ode/tests/test_radau5_van_der_pol.rs
+++ b/russell_ode/tests/test_radau5_van_der_pol.rs
@@ -16,7 +16,7 @@ fn test_radau5_van_der_pol() {
     out.enable_dense(0.2, &[0, 1]).unwrap();
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver
         .solve(&mut data.y0, data.x0, data.x1, None, Some(&mut out), &mut args)
         .unwrap();
diff --git a/russell_ode/tests/test_radau5_van_der_pol_debug.rs b/russell_ode/tests/test_radau5_van_der_pol_debug.rs
index 9f3b3448..fd7058c4 100644
--- a/russell_ode/tests/test_radau5_van_der_pol_debug.rs
+++ b/russell_ode/tests/test_radau5_van_der_pol_debug.rs
@@ -13,7 +13,7 @@ fn test_radau5_van_der_pol_debug() {
     params.debug = true;
 
     // solve the ODE system
-    let mut solver = OdeSolver::new(params, system).unwrap();
+    let mut solver = OdeSolver::new(params, &system).unwrap();
     solver
         .solve(&mut data.y0, data.x0, data.x1, None, None, &mut args)
         .unwrap();