diff --git a/.gitignore b/.gitignore
index f10bc18e..81b620f0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -12,3 +12,4 @@ good.mtx
 bad.mtx
 good.smat
 bad.smat
+rustc-ice-*.txt
diff --git a/russell_ode/src/erk_dense_out.rs b/russell_ode/src/erk_dense_out.rs
index fb8682a9..6bbe1ca0 100644
--- a/russell_ode/src/erk_dense_out.rs
+++ b/russell_ode/src/erk_dense_out.rs
@@ -57,19 +57,16 @@ impl ErkDenseOut {
     }
 
     /// Updates the data and returns the number of function evaluations
-    pub(crate) fn update<'a, F, A>(
+    pub(crate) fn update<'a, A>(
         &mut self,
-        system: &System<F, A>,
+        system: &System<A>,
         x: f64,
         y: &Vector,
         h: f64,
         w: &Vector,
         k: &Vec<Vector>,
         args: &mut A,
-    ) -> Result<usize, StrError>
-    where
-        F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-    {
+    ) -> Result<usize, StrError> {
         let mut n_function_eval = 0;
 
         if self.method == Method::DoPri5 {
diff --git a/russell_ode/src/euler_backward.rs b/russell_ode/src/euler_backward.rs
index 544e05b7..b7067ac3 100644
--- a/russell_ode/src/euler_backward.rs
+++ b/russell_ode/src/euler_backward.rs
@@ -4,15 +4,12 @@ use russell_lab::{vec_copy, vec_rms_scaled, vec_update, Vector};
 use russell_sparse::{numerical_jacobian, LinSolver, SparseMatrix};
 
 /// Implements the backward Euler (implicit) solver (implicit, order 1, unconditionally stable)
-pub(crate) struct EulerBackward<'a, F, A>
-where
-    F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-{
+pub(crate) struct EulerBackward<'a, A> {
     /// Holds the parameters
     params: Params,
 
     /// ODE system
-    system: &'a System<'a, F, A>,
+    system: &'a System<'a, A>,
 
     /// Vector holding the function evaluation
     ///
@@ -35,12 +32,9 @@ where
     solver: LinSolver<'a>,
 }
 
-impl<'a, F, A> EulerBackward<'a, F, A>
-where
-    F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-{
+impl<'a, A> EulerBackward<'a, A> {
     /// Allocates a new instance
-    pub fn new(params: Params, system: &'a System<'a, F, A>) -> Self {
+    pub fn new(params: Params, system: &'a System<'a, A>) -> Self {
         let ndim = system.ndim;
         let jac_nnz = if params.newton.use_numerical_jacobian {
             ndim * ndim
@@ -61,10 +55,7 @@ where
     }
 }
 
-impl<'a, F, A> OdeSolverTrait<A> for EulerBackward<'a, F, A>
-where
-    F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-{
+impl<'a, A> OdeSolverTrait<A> for EulerBackward<'a, A> {
     /// Enables dense output
     fn enable_dense_output(&mut self) -> Result<(), StrError> {
         Err("dense output is not available for the BwEuler method")
diff --git a/russell_ode/src/euler_forward.rs b/russell_ode/src/euler_forward.rs
index 42b8f3c9..dc8ba00d 100644
--- a/russell_ode/src/euler_forward.rs
+++ b/russell_ode/src/euler_forward.rs
@@ -6,12 +6,9 @@ use russell_lab::{vec_add, vec_copy, Vector};
 ///
 /// **Warning:** This method is interesting for didactic purposes only
 /// and should not be used in production codes.
-pub(crate) struct EulerForward<'a, F, A>
-where
-    F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-{
+pub(crate) struct EulerForward<'a, A> {
     /// ODE system
-    system: &'a System<'a, F, A>,
+    system: &'a System<'a, A>,
 
     /// Vector holding the function evaluation
     ///
@@ -22,12 +19,9 @@ where
     w: Vector,
 }
 
-impl<'a, F, A> EulerForward<'a, F, A>
-where
-    F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-{
+impl<'a, A> EulerForward<'a, A> {
     /// Allocates a new instance
-    pub fn new(system: &'a System<'a, F, A>) -> Self {
+    pub fn new(system: &'a System<'a, A>) -> Self {
         let ndim = system.ndim;
         EulerForward {
             system,
@@ -37,10 +31,7 @@ where
     }
 }
 
-impl<'a, F, A> OdeSolverTrait<A> for EulerForward<'a, F, A>
-where
-    F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-{
+impl<'a, A> OdeSolverTrait<A> for EulerForward<'a, A> {
     /// Enables dense output
     fn enable_dense_output(&mut self) -> Result<(), StrError> {
         Err("dense output is not available for the FwEuler method")
diff --git a/russell_ode/src/explicit_runge_kutta.rs b/russell_ode/src/explicit_runge_kutta.rs
index b0829fdb..061e9885 100644
--- a/russell_ode/src/explicit_runge_kutta.rs
+++ b/russell_ode/src/explicit_runge_kutta.rs
@@ -20,15 +20,12 @@ use russell_lab::{format_fortran, vec_copy, vec_update, Matrix, Vector};
 /// 2. E. Hairer, G. Wanner (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(crate) struct ExplicitRungeKutta<'a, F, A>
-where
-    F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-{
+pub(crate) struct ExplicitRungeKutta<'a, A> {
     /// Holds the parameters
     params: Params,
 
     /// ODE system
-    system: &'a System<'a, F, A>,
+    system: &'a System<'a, A>,
 
     /// Information such as implicit, embedded, etc.
     info: Information,
@@ -78,12 +75,9 @@ where
     dense_out: Option<ErkDenseOut>,
 }
 
-impl<'a, F, A> ExplicitRungeKutta<'a, F, A>
-where
-    F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-{
+impl<'a, A> ExplicitRungeKutta<'a, A> {
     /// Allocates a new instance
-    pub fn new(params: Params, system: &'a System<'a, F, A>) -> Result<Self, StrError> {
+    pub fn new(params: Params, system: &'a System<'a, A>) -> Result<Self, StrError> {
         // Runge-Kutta coefficients
         #[rustfmt::skip]
         let (aa, bb, cc) = match params.method {
@@ -153,10 +147,7 @@ where
     }
 }
 
-impl<'a, F, A> OdeSolverTrait<A> for ExplicitRungeKutta<'a, F, A>
-where
-    F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-{
+impl<'a, A> OdeSolverTrait<A> for ExplicitRungeKutta<'a, A> {
     /// Enables dense output
     fn enable_dense_output(&mut self) -> Result<(), StrError> {
         self.dense_out = Some(ErkDenseOut::new(self.params.method, self.system.ndim)?);
diff --git a/russell_ode/src/ode_solver.rs b/russell_ode/src/ode_solver.rs
index 67c6d5f0..0275a0f7 100644
--- a/russell_ode/src/ode_solver.rs
+++ b/russell_ode/src/ode_solver.rs
@@ -133,13 +133,10 @@ impl<'a, A> OdeSolver<'a, A> {
     ///
     /// # Generics
     ///
-    /// The generic arguments here are:
-    ///
-    /// * `F` -- function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)`
-    /// * `A` -- generic argument to assist in the `F` and Jacobian functions. It may be simply [crate::NoArgs] indicating that no arguments are needed.
-    pub fn new<F>(params: Params, system: &'a System<'a, F, A>) -> Result<Self, StrError>
+    /// * `A` -- generic argument to assist in the f(x,y) and Jacobian functions.
+    ///   It may be simply [NoArgs] indicating that no arguments are needed.
+    pub fn new(params: Params, system: &'a System<'a, A>) -> Result<Self, StrError>
     where
-        F: 'a + Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
         A: 'a,
     {
         if system.mass_matrix.is_some() && params.method != Method::Radau5 {
diff --git a/russell_ode/src/output.rs b/russell_ode/src/output.rs
index 8cbdf71f..22a1eb14 100644
--- a/russell_ode/src/output.rs
+++ b/russell_ode/src/output.rs
@@ -38,10 +38,8 @@ pub struct OutCount {
 ///
 /// # Generics
 ///
-/// The generic arguments are:
-///
-/// * `A` -- Is auxiliary argument for the `F`, `J`, `YxFunction`, and `OutCallback` functions.
-///   It may be simply [crate::NoArgs] indicating that no arguments are effectively used.
+/// * `A` -- generic argument to assist in the f(x,y) and Jacobian functions.
+///   It may be simply [crate::NoArgs] indicating that no arguments are needed.
 pub struct Output<'a, A> {
     /// Indicates whether the solver called initialize or not
     initialized: bool,
diff --git a/russell_ode/src/radau5.rs b/russell_ode/src/radau5.rs
index ab09c1a0..54dbda1c 100644
--- a/russell_ode/src/radau5.rs
+++ b/russell_ode/src/radau5.rs
@@ -24,15 +24,12 @@ use std::thread;
 /// 2. E. Hairer, G. Wanner (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(crate) struct Radau5<'a, F, A>
-where
-    F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-{
+pub(crate) struct Radau5<'a, A> {
     /// Holds the parameters
     params: Params,
 
     /// ODE system
-    system: &'a System<'a, F, A>,
+    system: &'a System<'a, A>,
 
     /// Holds the Jacobian matrix. J = df/dy
     jj: SparseMatrix,
@@ -118,12 +115,9 @@ where
     dw12: ComplexVector, // packed (dw1, dw2)
 }
 
-impl<'a, F, A> Radau5<'a, F, A>
-where
-    F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-{
+impl<'a, A> Radau5<'a, A> {
     /// Allocates a new instance
-    pub fn new(params: Params, system: &'a System<'a, F, A>) -> Self {
+    pub fn new(params: Params, system: &'a System<'a, A>) -> Self {
         let ndim = system.ndim;
         let mass_nnz = match system.mass_matrix.as_ref() {
             Some(mass) => mass.get_info().2,
@@ -331,10 +325,7 @@ where
     }
 }
 
-impl<'a, F, A> OdeSolverTrait<A> for Radau5<'a, F, A>
-where
-    F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-{
+impl<'a, A> OdeSolverTrait<A> for Radau5<'a, A> {
     /// Enables dense output
     fn enable_dense_output(&mut self) -> Result<(), StrError> {
         Ok(())
diff --git a/russell_ode/src/samples.rs b/russell_ode/src/samples.rs
index 36ed3de8..ca8b7da5 100644
--- a/russell_ode/src/samples.rs
+++ b/russell_ode/src/samples.rs
@@ -1,4 +1,3 @@
-use crate::StrError;
 use crate::{NoArgs, PdeDiscreteLaplacian2d, Side, System};
 use russell_lab::math::PI;
 use russell_lab::Vector;
@@ -42,7 +41,7 @@ impl Samples {
     /// * `args` -- is a placeholder variable with the arguments to F and J
     /// * `y_fn_x` -- is a function to compute the analytical solution
     pub fn simple_equation_constant<'a>() -> (
-        System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
+        System<'a, NoArgs>,
         f64,
         Vector,
         NoArgs,
@@ -153,7 +152,7 @@ impl Samples {
         symmetric: bool,
         genie: Genie,
     ) -> (
-        System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
+        System<'a, NoArgs>,
         f64,
         Vector,
         NoArgs,
@@ -257,13 +256,7 @@ impl Samples {
     /// * Hairer E, Nørsett, SP, Wanner G (2008) Solving Ordinary Differential Equations I.
     ///   Non-stiff Problems. Second Revised Edition. Corrected 3rd printing 2008. Springer Series
     ///   in Computational Mathematics, 528p
-    pub fn brusselator_ode<'a>() -> (
-        System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
-        f64,
-        Vector,
-        NoArgs,
-        Vector,
-    ) {
+    pub fn brusselator_ode<'a>() -> (System<'a, NoArgs>, f64, Vector, NoArgs, Vector) {
         // system
         let ndim = 2;
         let jac_nnz = 4;
@@ -502,16 +495,7 @@ impl Samples {
         npoint: usize,
         second_book: bool,
         ignore_diffusion: bool,
-    ) -> (
-        System<
-            'a,
-            impl Fn(&mut Vector, f64, &Vector, &mut PdeDiscreteLaplacian2d) -> Result<(), StrError>,
-            PdeDiscreteLaplacian2d,
-        >,
-        f64,
-        Vector,
-        PdeDiscreteLaplacian2d,
-    ) {
+    ) -> (System<'a, PdeDiscreteLaplacian2d>, f64, Vector, PdeDiscreteLaplacian2d) {
         // constants
         let (kx, ky) = (alpha, alpha);
         let (xmin, xmax) = (0.0, 1.0);
@@ -663,14 +647,7 @@ impl Samples {
     /// * Hairer E, Nørsett, SP, Wanner G (2008) Solving Ordinary Differential Equations I.
     ///   Non-stiff Problems. Second Revised Edition. Corrected 3rd printing 2008. Springer Series
     ///   in Computational Mathematics, 528p
-    pub fn arenstorf<'a>() -> (
-        System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
-        f64,
-        Vector,
-        f64,
-        NoArgs,
-        Vector,
-    ) {
+    pub fn arenstorf<'a>() -> (System<'a, NoArgs>, f64, Vector, f64, NoArgs, Vector) {
         // constants
         const MU: f64 = 0.012277471;
         const MD: f64 = 1.0 - MU;
@@ -788,7 +765,7 @@ impl Samples {
     ///   Stiff and Differential-Algebraic Problems. Second Revised Edition.
     ///   Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p
     pub fn hairer_wanner_eq1<'a>() -> (
-        System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
+        System<'a, NoArgs>,
         f64,
         Vector,
         NoArgs,
@@ -861,12 +838,7 @@ 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>() -> (
-        System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
-        f64,
-        Vector,
-        NoArgs,
-    ) {
+    pub fn robertson<'a>() -> (System<'a, NoArgs>, f64, Vector, NoArgs) {
         // system
         let ndim = 3;
         let mut system = System::new(ndim, |f: &mut Vector, _x: f64, y: &Vector, _args: &mut NoArgs| {
@@ -942,16 +914,7 @@ 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>(
-        epsilon: f64,
-        stationary: bool,
-    ) -> (
-        System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
-        f64,
-        Vector,
-        f64,
-        NoArgs,
-    ) {
+    pub fn van_der_pol<'a>(epsilon: f64, stationary: bool) -> (System<'a, NoArgs>, f64, Vector, f64, NoArgs) {
         // constants
         let x0 = 0.0;
         let mut y0 = Vector::from(&[2.0, -0.6]);
@@ -1071,12 +1034,7 @@ 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 amplifier1t<'a>() -> (
-        System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
-        f64,
-        Vector,
-        NoArgs,
-    ) {
+    pub fn amplifier1t<'a>() -> (System<'a, NoArgs>, f64, Vector, NoArgs) {
         // constants
         const ALPHA: f64 = 0.99;
         const GAMMA: f64 = 1.0 - ALPHA;
@@ -1174,7 +1132,7 @@ 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>() -> (
-        System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
+        System<'a, NoArgs>,
         f64,
         Vector,
         NoArgs,
@@ -1250,7 +1208,7 @@ 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>() -> (
-        System<'a, impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, NoArgs>,
+        System<'a, NoArgs>,
         f64,
         Vector,
         NoArgs,
diff --git a/russell_ode/src/system.rs b/russell_ode/src/system.rs
index ab03a088..14354ed0 100644
--- a/russell_ode/src/system.rs
+++ b/russell_ode/src/system.rs
@@ -37,10 +37,8 @@ pub type NoArgs = u8;
 ///
 /// # Generics
 ///
-/// The generic arguments here are:
-///
-/// * `F` -- function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)`
-/// * `A` -- generic argument to assist in the `F` and Jacobian functions. It may be simply [NoArgs] indicating that no arguments are needed.
+/// * `A` -- generic argument to assist in the f(x,y) and Jacobian functions.
+///   It may be simply [NoArgs] indicating that no arguments are needed.
 ///
 /// # Important
 ///
@@ -63,15 +61,12 @@ pub type NoArgs = u8;
 /// 2. E. Hairer, G. Wanner (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 struct System<'a, F, A>
-where
-    F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-{
+pub struct System<'a, A> {
     /// System dimension
     pub(crate) ndim: usize,
 
     /// ODE system function
-    pub(crate) function: F,
+    pub(crate) function: Box<dyn Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError> + 'a>,
 
     /// Jacobian function
     pub(crate) jacobian: Option<Box<dyn Fn(&mut CooMatrix, f64, f64, &Vector, &mut A) -> Result<(), StrError> + 'a>>,
@@ -95,10 +90,7 @@ where
     phantom: PhantomData<fn() -> A>,
 }
 
-impl<'a, F, A> System<'a, F, A>
-where
-    F: Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError>,
-{
+impl<'a, A> System<'a, A> {
     /// Allocates a new instance
     ///
     /// # Input
@@ -111,10 +103,8 @@ where
     ///
     /// # Generics
     ///
-    /// The generic arguments here are:
-    ///
-    /// * `F` -- function to compute the `f` vector: `(f: &mut Vector, x: f64, y: &Vector, args: &mut A)`
-    /// * `A` -- generic argument to assist in the `F` and Jacobian functions. It may be simply [NoArgs] indicating that no arguments are needed.
+    /// * `A` -- generic argument to assist in the f(x,y) and Jacobian functions.
+    ///   It may be simply [NoArgs] indicating that no arguments are needed.
     ///
     /// # Examples
     ///
@@ -161,10 +151,10 @@ where
     ///     Ok(())
     /// }
     /// ```
-    pub fn new(ndim: usize, function: F) -> Self {
+    pub fn new(ndim: usize, function: impl Fn(&mut Vector, f64, &Vector, &mut A) -> Result<(), StrError> + 'a) -> Self {
         System {
             ndim,
-            function,
+            function: Box::new(function),
             jacobian: None,
             jac_nnz: ndim * ndim,
             sym_jac: None,