diff --git a/Cargo.toml b/Cargo.toml
index 4df84cfbc..4623b7acb 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -28,7 +28,7 @@ tracing-subscriber = { version = "0.3.19", features = [
     "time",
 ] }
 faer = "0.24.0"
-pharmsol = "=0.24.1"
+pharmsol = "=0.25"
 rand = "0.9.0"
 anyhow = "1.0.100"
 rayon = "1.10.0"
diff --git a/benches/bimodal_ke.rs b/benches/bimodal_ke.rs
index b8ce6c55e..fc0d5d01d 100644
--- a/benches/bimodal_ke.rs
+++ b/benches/bimodal_ke.rs
@@ -18,6 +18,7 @@ fn create_equation() -> equation::ODE {
             y[1] = x[0] / v;
         },
     )
+    .with_default_cache()
 }
 
 fn create_parameters() -> Parameters {
diff --git a/examples/bestdose.rs b/examples/bestdose.rs
index 6968514a6..3fa34c8fc 100644
--- a/examples/bestdose.rs
+++ b/examples/bestdose.rs
@@ -16,7 +16,8 @@ fn main() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    };
+    }
+    .with_default_cache();
 
     let params = Parameters::new()
         .add("ke", 0.001, 3.0)
diff --git a/examples/bestdose_auc.rs b/examples/bestdose_auc.rs
index 70aecb429..475b72562 100644
--- a/examples/bestdose_auc.rs
+++ b/examples/bestdose_auc.rs
@@ -19,7 +19,8 @@ fn main() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    };
+    }
+    .with_default_cache();
 
     // Minimal parameter ranges
     let params = Parameters::new()
diff --git a/examples/bestdose_bounds.rs b/examples/bestdose_bounds.rs
index 4b2eab3bb..7c57fbd4a 100644
--- a/examples/bestdose_bounds.rs
+++ b/examples/bestdose_bounds.rs
@@ -19,7 +19,8 @@ fn main() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    };
+    }
+    .with_default_cache();
 
     let params = Parameters::new()
         .add("ke", 0.001, 3.0)
diff --git a/examples/bimodal_ke/main.rs b/examples/bimodal_ke/main.rs
index 0ea9bb68c..58e2d008b 100644
--- a/examples/bimodal_ke/main.rs
+++ b/examples/bimodal_ke/main.rs
@@ -13,18 +13,17 @@ fn main() -> Result<()> {
             y[1] = x[0] / v;
         },
     }
-    .with_solver(OdeSolver::ExplicitRk(ExplicitRkTableau::Tsit45));
+    .with_solver(OdeSolver::ExplicitRk(ExplicitRkTableau::Tsit45))
+    .with_default_cache();
 
     let params = Parameters::new()
         .add("ke", 0.001, 3.0)
         .add("v", 25.0, 250.0);
 
-    let ems = AssayErrorModels::new()
-        .add(
-            1,
-            AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
-        )
-        .unwrap();
+    let ems = AssayErrorModels::new().add(
+        1,
+        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.5, 0.0, 0.0), 0.0),
+    )?;
 
     let mut settings = Settings::builder()
         .set_algorithm(Algorithm::NPAG)
diff --git a/examples/drusano/main.rs b/examples/drusano/main.rs
index 640fb20eb..58ec866ed 100644
--- a/examples/drusano/main.rs
+++ b/examples/drusano/main.rs
@@ -69,7 +69,8 @@ fn main() -> Result<()> {
             y[3] = x[3].log10();
             y[4] = x[4].log10();
         },
-    };
+    }
+    .with_default_cache();
 
     let params = Parameters::new()
         .add("v1", 5.0, 160.0)
diff --git a/examples/iov/main.rs b/examples/iov/main.rs
index 222015dfc..53abf12db 100644
--- a/examples/iov/main.rs
+++ b/examples/iov/main.rs
@@ -29,7 +29,8 @@ fn main() -> Result<()> {
     )
     .with_nstates(2)
     .with_ndrugs(1)
-    .with_nout(1);
+    .with_nout(1)
+    .with_default_cache();
 
     let params = Parameters::new().add("ke0", 0.001, 2.0);
 
diff --git a/examples/meta/main.rs b/examples/meta/main.rs
index 61dc77ebf..46f7fed41 100644
--- a/examples/meta/main.rs
+++ b/examples/meta/main.rs
@@ -26,7 +26,8 @@ fn main() {
             y[1] = x[0] / v;
             y[2] = x[1] / v2;
         },
-    };
+    }
+    .with_default_cache();
 
     let params = Parameters::new()
         .add("cls", 0.1, 10.0)
diff --git a/examples/neely/main.rs b/examples/neely/main.rs
index 3cf903b19..a61f700ca 100644
--- a/examples/neely/main.rs
+++ b/examples/neely/main.rs
@@ -45,7 +45,9 @@ fn main() {
             y[2] = x[2] / vm1;
             y[3] = x[3] / vm2;
         },
-    };
+    }
+    .with_default_cache();
+
     let params = Parameters::new()
         .add("cls", 0.0, 0.4)
         .add("k30", 0.0, 0.5)
diff --git a/examples/new_iov/main.rs b/examples/new_iov/main.rs
index 40490b7f2..d50a9a970 100644
--- a/examples/new_iov/main.rs
+++ b/examples/new_iov/main.rs
@@ -29,7 +29,8 @@ fn main() {
     )
     .with_nstates(2)
     .with_ndrugs(1)
-    .with_nout(1);
+    .with_nout(1)
+    .with_default_cache();
 
     let params = Parameters::new()
         .add("ke0", 0.0001, 2.4)
@@ -49,7 +50,6 @@ fn main() {
         .build();
 
     settings.set_cycles(1000);
-    settings.set_cache(true);
     settings.set_output_path("examples/new_iov/output");
     settings.set_prior(Prior::sobol(100, 347));
 
diff --git a/examples/theophylline/main.rs b/examples/theophylline/main.rs
index 9addb5e2f..8eb78e3b4 100644
--- a/examples/theophylline/main.rs
+++ b/examples/theophylline/main.rs
@@ -11,7 +11,8 @@ fn main() {
             fetch_params!(p, _ka, _ke, v);
             y[0] = x[1] * 1000.0 / v;
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new()
         .add("ka", 0.001, 3.0)
diff --git a/examples/two_eq_lag/main.rs b/examples/two_eq_lag/main.rs
index 92454a6b6..ed4dd2fc1 100644
--- a/examples/two_eq_lag/main.rs
+++ b/examples/two_eq_lag/main.rs
@@ -20,7 +20,8 @@ fn main() {
             fetch_params!(p, _ka, _ke, _tlag, v);
             y[1] = x[1] / v;
         },
-    };
+    }
+    .with_default_cache();
     // let eq = Equation::new_analytical(
     //     one_compartment_with_absorption,
     //     |_p, _cov| {},
diff --git a/examples/vanco.rs b/examples/vanco.rs
index f76b60c67..7e111f34a 100644
--- a/examples/vanco.rs
+++ b/examples/vanco.rs
@@ -12,7 +12,8 @@ fn main() {
         out: |x, _p, _t, _cov, y| {
             y[0] = x[1];
         },
-    };
+    }
+    .with_default_cache();
     // same eq but analytical
     // let eq = Equation::new_analytical(
     //     two_compartments,
diff --git a/examples/vanco_sde/main.rs b/examples/vanco_sde/main.rs
index 024659c4e..f600251d5 100644
--- a/examples/vanco_sde/main.rs
+++ b/examples/vanco_sde/main.rs
@@ -27,7 +27,8 @@ fn main() {
             y[0] = x[1] / (vol * wt);
         },
         100,
-    );
+    )
+    .with_default_cache();
 
     // let ode = equation::ODE::new(
     //     |x, p, _t, dx, _rateiv, _cov| {
@@ -69,7 +70,6 @@ fn main() {
         .build();
 
     settings.set_cycles(usize::MAX);
-    settings.set_cache(true);
     settings.set_output_path("examples/vanco_sde/output");
     settings.set_prior(Prior::sobol(100, 347));
     settings.initialize_logs().unwrap();
diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs
index 949ded9c0..0ce1dea8a 100644
--- a/src/algorithms/mod.rs
+++ b/src/algorithms/mod.rs
@@ -185,7 +185,7 @@ pub trait Algorithms<E: Equation + Send + 'static>: Sync + Send + 'static {
                     let outeqs = preds.iter().map(|x| x.outeq()).collect::<Vec<usize>>();
                     let states = preds
                         .iter()
-                        .map(|x| x.state().clone())
+                        .map(|x| x.state().to_vec())
                         .collect::<Vec<Vec<f64>>>();
 
                     tracing::debug!("\t\tTimes: {:?}", times);
diff --git a/src/routines/settings.rs b/src/routines/settings.rs
index 45383b287..e03fa1e31 100644
--- a/src/routines/settings.rs
+++ b/src/routines/settings.rs
@@ -84,10 +84,6 @@ impl Settings {
         self.config.algorithm = algorithm;
     }
 
-    pub fn set_cache(&mut self, cache: bool) {
-        self.config.cache = cache;
-    }
-
     pub fn set_idelta(&mut self, idelta: f64) {
         self.predictions.idelta = idelta;
     }
@@ -148,8 +144,6 @@ pub struct Config {
     pub cycles: usize,
     /// Denotes the algorithm to use
     pub algorithm: Algorithm,
-    /// If true (default), cache predicted values
-    pub cache: bool,
     /// Should a progress bar be displayed for the first cycle
     ///
     /// The progress bar is not written to logs, but is written to stdout. It incurs a minor performance penalty.
@@ -161,7 +155,6 @@ impl Default for Config {
         Config {
             cycles: 100,
             algorithm: Algorithm::NPAG,
-            cache: true,
             progress: true,
         }
     }
@@ -603,7 +596,6 @@ mod tests {
 
         assert_eq!(settings.config.algorithm, Algorithm::NPAG);
         assert_eq!(settings.config.cycles, 100);
-        assert_eq!(settings.config.cache, true);
         assert_eq!(settings.parameters().names(), vec!["Ke", "V"]);
     }
 }
diff --git a/tests/bestdose_tests.rs b/tests/bestdose_tests.rs
index 145bf11ae..57118c86a 100644
--- a/tests/bestdose_tests.rs
+++ b/tests/bestdose_tests.rs
@@ -21,7 +21,8 @@ fn test_infusion_mask_inclusion() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
 
@@ -136,7 +137,8 @@ fn test_fixed_infusion_preservation() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new()
         .add("ke", 0.001, 3.0)
@@ -229,7 +231,8 @@ fn test_dose_count_validation() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
     let ems = AssayErrorModels::new().add(
@@ -299,7 +302,8 @@ fn test_empty_observations_validation() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
     let ems = AssayErrorModels::new().add(
@@ -359,7 +363,8 @@ fn test_basic_auc_mode() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
 
@@ -444,7 +449,8 @@ fn test_infusion_auc_mode() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
 
@@ -548,7 +554,8 @@ fn test_multi_outeq_auc_mode() -> Result<()> {
             y[0] = x[0] / v; // outeq 0: concentration
             y[1] = x[0]; // outeq 1: amount
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
 
@@ -618,7 +625,8 @@ fn test_multi_outeq_auc_optimization() -> Result<()> {
             y[0] = x[0] / v;
             y[1] = x[0];
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
     let error_model = AssayErrorModel::additive(ErrorPoly::new(0.0, 5.0, 0.0, 0.0), 0.0);
@@ -703,7 +711,8 @@ fn test_auc_from_zero_single_dose() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new().add("ke", 0.2, 0.4).add("v", 40.0, 60.0);
 
@@ -787,7 +796,8 @@ fn test_auc_from_last_dose_maintenance() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new().add("ke", 0.2, 0.4).add("v", 40.0, 60.0);
 
@@ -875,7 +885,8 @@ fn test_auc_modes_comparison() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new().add("ke", 0.3, 0.3).add("v", 50.0, 50.0);
 
@@ -1004,7 +1015,8 @@ fn test_auc_from_last_dose_multiple_observations() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new().add("ke", 0.2, 0.4).add("v", 40.0, 60.0);
 
@@ -1098,7 +1110,8 @@ fn test_auc_from_last_dose_no_prior_dose() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new().add("ke", 0.2, 0.4).add("v", 40.0, 60.0);
 
@@ -1187,7 +1200,8 @@ fn test_dose_range_bounds_respected() -> Result<()> {
             fetch_params!(p, _ke, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     let params = Parameters::new().add("ke", 0.1, 0.5).add("v", 40.0, 60.0);
 
@@ -1281,6 +1295,7 @@ fn one_compartment_model() -> pharmsol::ODE {
             y[0] = x[0] / v;
         },
     )
+    .with_default_cache()
 }
 
 /// Helper to build minimal settings for tests (no posterior refinement)
diff --git a/tests/onecomp.rs b/tests/onecomp.rs
index 58230aca0..3f20ed9b4 100644
--- a/tests/onecomp.rs
+++ b/tests/onecomp.rs
@@ -16,7 +16,8 @@ fn test_one_compartment_npag() -> Result<()> {
             fetch_params!(p, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     // Define parameters
     let params = Parameters::new().add("ke", 0.1, 1.0).add("v", 1.0, 20.0);
@@ -83,7 +84,8 @@ fn test_one_compartment_npod() -> Result<()> {
             fetch_params!(p, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     // Define parameters
     let params = Parameters::new().add("ke", 0.1, 1.0).add("v", 1.0, 20.0);
@@ -150,7 +152,8 @@ fn test_one_compartment_postprob() -> Result<()> {
             fetch_params!(p, v);
             y[0] = x[0] / v;
         },
-    );
+    )
+    .with_default_cache();
 
     // Define parameters
     let params = Parameters::new().add("ke", 0.1, 1.0).add("v", 1.0, 20.0);
diff --git a/tests/settings_tests.rs b/tests/settings_tests.rs
index 23992c34b..bcb848b4a 100644
--- a/tests/settings_tests.rs
+++ b/tests/settings_tests.rs
@@ -97,10 +97,6 @@ fn test_settings_setters() -> Result<()> {
     settings.set_algorithm(Algorithm::NPOD);
     assert_eq!(settings.config().algorithm, Algorithm::NPOD);
 
-    // Test set_cache
-    settings.set_cache(false);
-    assert_eq!(settings.config().cache, false);
-
     // Test set_idelta
     settings.set_idelta(0.5);
     assert_eq!(settings.predictions().idelta, 0.5);
@@ -212,7 +208,6 @@ fn test_config_accessors() -> Result<()> {
     // Test default values
     assert_eq!(config.algorithm, Algorithm::NPAG);
     assert!(config.cycles > 0);
-    assert_eq!(config.cache, true);
 
     Ok(())
 }