diff --git a/src/lib.rs b/src/lib.rs
index 981e1f00..c12599e5 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -137,7 +137,7 @@ pub mod solver;
 pub mod vector;
 
 use linear_solver::LinearSolver;
-pub use linear_solver::NalgebraLU;
+pub use linear_solver::{FaerLU, NalgebraLU};
 
 #[cfg(feature = "sundials")]
 pub use matrix::sundials::SundialsMatrix;
@@ -181,7 +181,7 @@ mod tests {
             .p([0.04, 1.0e4, 3.0e7])
             .rtol(1e-4)
             .atol([1.0e-8, 1.0e-6, 1.0e-6])
-            .build_ode(
+            .build_ode_dense(
                 |x: &V, p: &V, _t: T, y: &mut V| {
                     y[0] = -p[0] * x[0] + p[1] * x[1] * x[2];
                     y[1] = p[0] * x[0] - p[1] * x[1] * x[2] - p[2] * x[1] * x[1];
@@ -217,11 +217,12 @@ mod tests {
     fn test_readme_faer() {
         type T = f64;
         type V = faer::Col<f64>;
+        type M = faer::Mat<f64>;
         let problem = OdeBuilder::new()
             .p([0.04, 1.0e4, 3.0e7])
             .rtol(1e-4)
             .atol([1.0e-8, 1.0e-6, 1.0e-6])
-            .build_ode(
+            .build_ode_dense(
                 |x: &V, p: &V, _t: T, y: &mut V| {
                     y[0] = -p[0] * x[0] + p[1] * x[1] * x[2];
                     y[1] = p[0] * x[0] - p[1] * x[1] * x[2] - p[2] * x[1] * x[1];
@@ -239,7 +240,7 @@ mod tests {
             )
             .unwrap();
 
-        let mut solver = Bdf::default();
+        let mut solver = Bdf::<M, _>::default();
 
         let t = 0.4;
         let y = solver.solve(&problem, t).unwrap();
diff --git a/src/linear_solver/sundials.rs b/src/linear_solver/sundials.rs
index a964999d..175af9e9 100644
--- a/src/linear_solver/sundials.rs
+++ b/src/linear_solver/sundials.rs
@@ -23,6 +23,15 @@ where
     matrix: SundialsMatrix,
 }
 
+impl<Op> Default for SundialsLinearSolver<Op>
+where
+    Op: LinearOp<M = SundialsMatrix, V = SundialsVector, T = realtype>,
+{
+    fn default() -> Self {
+        Self::new_dense()
+    }
+}
+
 impl<Op> SundialsLinearSolver<Op>
 where
     Op: LinearOp<M = SundialsMatrix, V = SundialsVector, T = realtype>,
diff --git a/src/matrix/default_solver.rs b/src/matrix/default_solver.rs
new file mode 100644
index 00000000..3855ae2f
--- /dev/null
+++ b/src/matrix/default_solver.rs
@@ -0,0 +1,10 @@
+use crate::{linear_solver::LinearSolver, op::LinearOp};
+
+use super::Matrix;
+
+pub trait DefaultSolver: Matrix {
+    type LS<C: LinearOp<M = Self, V = Self::V, T = Self::T>>: LinearSolver<C> + Default;
+    fn default_solver<C: LinearOp<M = Self, V = Self::V, T = Self::T>>() -> Self::LS<C> {
+        Self::LS::default()
+    }
+}
diff --git a/src/matrix/dense_faer_serial.rs b/src/matrix/dense_faer_serial.rs
index a450da26..fa816988 100644
--- a/src/matrix/dense_faer_serial.rs
+++ b/src/matrix/dense_faer_serial.rs
@@ -1,10 +1,16 @@
 use std::ops::{Mul, MulAssign};
 
+use super::default_solver::DefaultSolver;
 use super::{DenseMatrix, Matrix, MatrixCommon, MatrixView, MatrixViewMut};
 use crate::scalar::{IndexType, Scalar, Scale};
+use crate::{op::LinearOp, FaerLU};
 use anyhow::Result;
 use faer::{linalg::matmul::matmul, Col, ColMut, ColRef, Mat, MatMut, MatRef, Parallelism};
 
+impl<T: Scalar> DefaultSolver for Mat<T> {
+    type LS<C: LinearOp<M = Mat<T>, V = Col<T>, T = T>> = FaerLU<T, C>;
+}
+
 macro_rules! impl_matrix_common {
     ($mat_type:ty) => {
         impl<'a, T: Scalar> MatrixCommon for $mat_type {
diff --git a/src/matrix/dense_nalgebra_serial.rs b/src/matrix/dense_nalgebra_serial.rs
index ddc881f2..029f37da 100644
--- a/src/matrix/dense_nalgebra_serial.rs
+++ b/src/matrix/dense_nalgebra_serial.rs
@@ -3,9 +3,16 @@ use std::ops::{Mul, MulAssign};
 use anyhow::Result;
 use nalgebra::{DMatrix, DMatrixView, DMatrixViewMut, DVector, DVectorView, DVectorViewMut};
 
+use crate::op::LinearOp;
 use crate::{scalar::Scale, IndexType, Scalar};
 
-use crate::{DenseMatrix, Matrix, MatrixCommon, MatrixView, MatrixViewMut};
+use crate::{DenseMatrix, Matrix, MatrixCommon, MatrixView, MatrixViewMut, NalgebraLU};
+
+use super::default_solver::DefaultSolver;
+
+impl<T: Scalar> DefaultSolver for DMatrix<T> {
+    type LS<C: LinearOp<M = DMatrix<T>, V = DVector<T>, T = T>> = NalgebraLU<T, C>;
+}
 
 macro_rules! impl_matrix_common {
     ($matrix_type:ty) => {
diff --git a/src/matrix/mod.rs b/src/matrix/mod.rs
index 280bb3cd..6f2fbc3a 100644
--- a/src/matrix/mod.rs
+++ b/src/matrix/mod.rs
@@ -12,7 +12,9 @@ mod dense_nalgebra_serial;
 #[cfg(feature = "faer")]
 mod dense_faer_serial;
 
+pub mod default_solver;
 mod sparse_serial;
+
 #[cfg(feature = "sundials")]
 pub mod sundials;
 
diff --git a/src/matrix/sundials.rs b/src/matrix/sundials.rs
index 54b8253a..e55fb3c3 100644
--- a/src/matrix/sundials.rs
+++ b/src/matrix/sundials.rs
@@ -12,12 +12,13 @@ use sundials_sys::{
 
 use crate::{
     ode_solver::sundials::sundials_check,
+    op::LinearOp,
     scalar::scale,
     vector::sundials::{get_suncontext, SundialsVector},
-    IndexType, Scale, Vector,
+    IndexType, Scale, SundialsLinearSolver, Vector,
 };
 
-use super::{Matrix, MatrixCommon};
+use super::{default_solver::DefaultSolver, Matrix, MatrixCommon};
 use anyhow::anyhow;
 
 #[derive(Debug)]
@@ -79,6 +80,11 @@ impl Display for SundialsMatrix {
     }
 }
 
+impl DefaultSolver for SundialsMatrix {
+    type LS<C: LinearOp<M = SundialsMatrix, V = SundialsVector, T = realtype>> =
+        SundialsLinearSolver<C>;
+}
+
 impl MatrixCommon for SundialsMatrix {
     type V = SundialsVector;
     type T = realtype;
diff --git a/src/ode_solver/bdf/faer.rs b/src/ode_solver/bdf/faer.rs
index 57beac46..83531a10 100644
--- a/src/ode_solver/bdf/faer.rs
+++ b/src/ode_solver/bdf/faer.rs
@@ -1,67 +1,3 @@
-use faer::{Col, Mat};
-
-use crate::{
-    linear_solver::FaerLU, op::bdf::BdfCallable, Bdf, NewtonNonlinearSolver, NonLinearSolver,
-    OdeEquations, Scalar, VectorRef,
-};
-
-impl<T: Scalar, Eqn: OdeEquations<T = T, V = Col<T>, M = Mat<T>> + 'static> Default
-    for Bdf<Mat<T>, Eqn>
-{
-    fn default() -> Self {
-        let n = 1;
-        let linear_solver = FaerLU::default();
-        let mut nonlinear_solver = Box::new(NewtonNonlinearSolver::<BdfCallable<Eqn>>::new(
-            linear_solver,
-        ));
-        nonlinear_solver.set_max_iter(Self::NEWTON_MAXITER);
-        Self {
-            ode_problem: None,
-            nonlinear_solver,
-            order: 1,
-            n_equal_steps: 0,
-            diff: Mat::zeros(n, Self::MAX_ORDER + 3),
-            diff_tmp: Mat::zeros(n, Self::MAX_ORDER + 3),
-            gamma: vec![T::from(1.0); Self::MAX_ORDER + 1],
-            alpha: vec![T::from(1.0); Self::MAX_ORDER + 1],
-            error_const: vec![T::from(1.0); Self::MAX_ORDER + 1],
-            u: Mat::zeros(Self::MAX_ORDER + 1, Self::MAX_ORDER + 1),
-            statistics: super::BdfStatistics::default(),
-            state: None,
-        }
-    }
-}
-
-// implement clone for bdf
-impl<T: Scalar, Eqn: OdeEquations<T = T, V = Col<T>, M = Mat<T>> + 'static> Clone
-    for Bdf<Mat<T>, Eqn>
-where
-    for<'b> &'b Col<T>: VectorRef<Col<T>>,
-{
-    fn clone(&self) -> Self {
-        let n = self.diff.nrows();
-        let linear_solver = FaerLU::default();
-        let mut nonlinear_solver = Box::new(NewtonNonlinearSolver::<BdfCallable<Eqn>>::new(
-            linear_solver,
-        ));
-        nonlinear_solver.set_max_iter(Self::NEWTON_MAXITER);
-        Self {
-            ode_problem: self.ode_problem.clone(),
-            nonlinear_solver,
-            order: self.order,
-            n_equal_steps: self.n_equal_steps,
-            diff: Mat::zeros(n, Self::MAX_ORDER + 3),
-            diff_tmp: Mat::zeros(n, Self::MAX_ORDER + 3),
-            gamma: self.gamma.clone(),
-            alpha: self.alpha.clone(),
-            error_const: self.error_const.clone(),
-            u: Mat::zeros(Self::MAX_ORDER + 1, Self::MAX_ORDER + 1),
-            statistics: self.statistics.clone(),
-            state: self.state.clone(),
-        }
-    }
-}
-
 #[cfg(test)]
 mod test {
     use crate::{
@@ -72,14 +8,14 @@ mod test {
     type M = faer::Mat<f64>;
     #[test]
     fn bdf_no_set_problem() {
-        test_no_set_problem::<M, _>(Bdf::default())
+        test_no_set_problem::<M, _>(Bdf::<M, _>::default())
     }
     #[test]
     fn bdf_take_state() {
-        test_take_state::<M, _>(Bdf::default())
+        test_take_state::<M, _>(Bdf::<M, _>::default())
     }
     #[test]
     fn bdf_test_interpolate() {
-        test_interpolate::<M, _>(Bdf::default())
+        test_interpolate::<M, _>(Bdf::<M, _>::default())
     }
 }
diff --git a/src/ode_solver/bdf/mod.rs b/src/ode_solver/bdf/mod.rs
index f41aac11..0e287fca 100644
--- a/src/ode_solver/bdf/mod.rs
+++ b/src/ode_solver/bdf/mod.rs
@@ -7,9 +7,13 @@ use num_traits::{One, Pow, Zero};
 use serde::Serialize;
 
 use crate::{
-    matrix::MatrixRef, op::bdf::BdfCallable, scalar::scale, DenseMatrix, IndexType, MatrixViewMut,
-    NonLinearSolver, OdeSolverMethod, OdeSolverProblem, OdeSolverState, Scalar, SolverProblem,
-    Vector, VectorRef, VectorView, VectorViewMut,
+    matrix::{default_solver::DefaultSolver, Matrix, MatrixRef},
+    op::bdf::BdfCallable,
+    scalar::scale,
+    vector::DefaultDenseMatrix,
+    DenseMatrix, IndexType, MatrixViewMut, NewtonNonlinearSolver, NonLinearSolver, OdeSolverMethod,
+    OdeSolverProblem, OdeSolverState, Scalar, SolverProblem, Vector, VectorRef, VectorView,
+    VectorViewMut,
 };
 
 pub mod faer;
@@ -69,32 +73,73 @@ pub struct Bdf<M: DenseMatrix<T = Eqn::T, V = Eqn::V>, Eqn: OdeEquations> {
     gamma: Vec<Eqn::T>,
     error_const: Vec<Eqn::T>,
     statistics: BdfStatistics<Eqn::T>,
-    state: Option<OdeSolverState<Eqn::M>>,
+    state: Option<OdeSolverState<Eqn::V>>,
 }
 
-// impl<Eqn: OdeEquations<T = f64, V = Col<f64>, M = faer::Mat<f64>>> Default for Bdf<Mat<f64>, Eqn> {
-//     fn default() -> Self {
-//         let n = 1;
-//         let linear_solver = LU::default();
-//         let mut nonlinear_solver = Box::new(NewtonNonlinearSolver::<BdfCallable<Eqn>>::new(
-//             linear_solver,
-//         ));
-//         nonlinear_solver.set_max_iter(Self::NEWTON_MAXITER);
-//         Self {
-//             ode_problem: None,
-//             nonlinear_solver,
-//             order: 1,
-//             n_equal_steps: 0,
-//             diff: Mat::zeros(n, Self::MAX_ORDER + 3), //DMatrix::<T>::zeros(n, Self::MAX_ORDER + 3),
-//             diff_tmp: Mat::zeros(n, Self::MAX_ORDER + 3),
-//             gamma: vec![f64::from(1.0); Self::MAX_ORDER + 1],
-//             alpha: vec![f64::from(1.0); Self::MAX_ORDER + 1],
-//             error_const: vec![f64::from(1.0); Self::MAX_ORDER + 1],
-//             u: Mat::zeros(Self::MAX_ORDER + 1, Self::MAX_ORDER + 1),
-//             statistics: BdfStatistics::default(),
-//         }
-//     }
-// }
+impl<Eqn> Default for Bdf<<Eqn::V as DefaultDenseMatrix>::M, Eqn>
+where
+    Eqn: OdeEquations + 'static,
+    Eqn::M: DefaultSolver,
+    Eqn::V: DefaultDenseMatrix,
+    for<'b> &'b Eqn::V: VectorRef<Eqn::V>,
+    for<'b> &'b Eqn::M: MatrixRef<Eqn::M>,
+{
+    fn default() -> Self {
+        let n = 1;
+        let linear_solver = Eqn::M::default_solver();
+        let mut nonlinear_solver = Box::new(NewtonNonlinearSolver::<BdfCallable<Eqn>>::new(
+            linear_solver,
+        ));
+        nonlinear_solver.set_max_iter(Self::NEWTON_MAXITER);
+        type M<V> = <V as DefaultDenseMatrix>::M;
+        Self {
+            ode_problem: None,
+            nonlinear_solver,
+            order: 1,
+            n_equal_steps: 0,
+            diff: <M<Eqn::V> as Matrix>::zeros(n, Self::MAX_ORDER + 3), //DMatrix::<T>::zeros(n, Self::MAX_ORDER + 3),
+            diff_tmp: <M<Eqn::V> as Matrix>::zeros(n, Self::MAX_ORDER + 3),
+            gamma: vec![Eqn::T::from(1.0); Self::MAX_ORDER + 1],
+            alpha: vec![Eqn::T::from(1.0); Self::MAX_ORDER + 1],
+            error_const: vec![Eqn::T::from(1.0); Self::MAX_ORDER + 1],
+            u: <M<Eqn::V> as Matrix>::zeros(Self::MAX_ORDER + 1, Self::MAX_ORDER + 1),
+            statistics: BdfStatistics::default(),
+            state: None,
+        }
+    }
+}
+
+impl<M, Eqn> Clone for Bdf<M, Eqn>
+where
+    M: DenseMatrix<T = Eqn::T, V = Eqn::V>,
+    Eqn: OdeEquations + 'static,
+    Eqn::M: DefaultSolver,
+    for<'b> &'b Eqn::V: VectorRef<Eqn::V>,
+    for<'b> &'b Eqn::M: MatrixRef<Eqn::M>,
+{
+    fn clone(&self) -> Self {
+        let n = self.diff.nrows();
+        let linear_solver = Eqn::M::default_solver();
+        let mut nonlinear_solver = Box::new(NewtonNonlinearSolver::<BdfCallable<Eqn>>::new(
+            linear_solver,
+        ));
+        nonlinear_solver.set_max_iter(Self::NEWTON_MAXITER);
+        Self {
+            ode_problem: self.ode_problem.clone(),
+            nonlinear_solver,
+            order: self.order,
+            n_equal_steps: self.n_equal_steps,
+            diff: M::zeros(n, Self::MAX_ORDER + 3),
+            diff_tmp: M::zeros(n, Self::MAX_ORDER + 3),
+            gamma: self.gamma.clone(),
+            alpha: self.alpha.clone(),
+            error_const: self.error_const.clone(),
+            u: M::zeros(Self::MAX_ORDER + 1, Self::MAX_ORDER + 1),
+            statistics: self.statistics.clone(),
+            state: self.state.clone(),
+        }
+    }
+}
 
 impl<M: DenseMatrix<T = Eqn::T, V = Eqn::V>, Eqn: OdeEquations> Bdf<M, Eqn>
 where
@@ -261,15 +306,15 @@ where
         self.ode_problem.as_ref()
     }
 
-    fn state(&self) -> Option<&OdeSolverState<Eqn::M>> {
+    fn state(&self) -> Option<&OdeSolverState<Eqn::V>> {
         self.state.as_ref()
     }
 
-    fn take_state(&mut self) -> Option<OdeSolverState<<Eqn>::M>> {
+    fn take_state(&mut self) -> Option<OdeSolverState<<Eqn>::V>> {
         Option::take(&mut self.state)
     }
 
-    fn set_problem(&mut self, state: OdeSolverState<Eqn::M>, problem: &OdeSolverProblem<Eqn>) {
+    fn set_problem(&mut self, state: OdeSolverState<Eqn::V>, problem: &OdeSolverProblem<Eqn>) {
         let mut state = state;
         self.ode_problem = Some(problem.clone());
         let nstates = problem.eqn.nstates();
diff --git a/src/ode_solver/bdf/nalgebra.rs b/src/ode_solver/bdf/nalgebra.rs
index 32647597..ed53a9b1 100644
--- a/src/ode_solver/bdf/nalgebra.rs
+++ b/src/ode_solver/bdf/nalgebra.rs
@@ -1,67 +1,3 @@
-use crate::{
-    linear_solver::NalgebraLU, op::bdf::BdfCallable, Bdf, NewtonNonlinearSolver, NonLinearSolver,
-    OdeEquations, Scalar, VectorRef,
-};
-
-use nalgebra::{DMatrix, DVector};
-
-impl<T: Scalar, Eqn: OdeEquations<T = T, V = DVector<T>, M = DMatrix<T>> + 'static> Default
-    for Bdf<DMatrix<T>, Eqn>
-{
-    fn default() -> Self {
-        let n = 1;
-        let linear_solver = NalgebraLU::default();
-        let mut nonlinear_solver = Box::new(NewtonNonlinearSolver::<BdfCallable<Eqn>>::new(
-            linear_solver,
-        ));
-        nonlinear_solver.set_max_iter(Self::NEWTON_MAXITER);
-        Self {
-            ode_problem: None,
-            nonlinear_solver,
-            order: 1,
-            n_equal_steps: 0,
-            diff: DMatrix::<T>::zeros(n, Self::MAX_ORDER + 3),
-            diff_tmp: DMatrix::<T>::zeros(n, Self::MAX_ORDER + 3),
-            gamma: vec![T::from(1.0); Self::MAX_ORDER + 1],
-            alpha: vec![T::from(1.0); Self::MAX_ORDER + 1],
-            error_const: vec![T::from(1.0); Self::MAX_ORDER + 1],
-            u: DMatrix::<T>::zeros(Self::MAX_ORDER + 1, Self::MAX_ORDER + 1),
-            statistics: super::BdfStatistics::default(),
-            state: None,
-        }
-    }
-}
-
-// implement clone for bdf
-impl<T: Scalar, Eqn: OdeEquations<T = T, V = DVector<T>, M = DMatrix<T>> + 'static> Clone
-    for Bdf<DMatrix<T>, Eqn>
-where
-    for<'b> &'b DVector<T>: VectorRef<DVector<T>>,
-{
-    fn clone(&self) -> Self {
-        let n = self.diff.nrows();
-        let linear_solver = NalgebraLU::default();
-        let mut nonlinear_solver = Box::new(NewtonNonlinearSolver::<BdfCallable<Eqn>>::new(
-            linear_solver,
-        ));
-        nonlinear_solver.set_max_iter(Self::NEWTON_MAXITER);
-        Self {
-            ode_problem: self.ode_problem.clone(),
-            nonlinear_solver,
-            order: self.order,
-            n_equal_steps: self.n_equal_steps,
-            diff: DMatrix::zeros(n, Self::MAX_ORDER + 3),
-            diff_tmp: DMatrix::zeros(n, Self::MAX_ORDER + 3),
-            gamma: self.gamma.clone(),
-            alpha: self.alpha.clone(),
-            error_const: self.error_const.clone(),
-            u: DMatrix::zeros(Self::MAX_ORDER + 1, Self::MAX_ORDER + 1),
-            statistics: self.statistics.clone(),
-            state: self.state.clone(),
-        }
-    }
-}
-
 #[cfg(test)]
 mod test {
     use crate::{
diff --git a/src/ode_solver/builder.rs b/src/ode_solver/builder.rs
index 89855260..775405d2 100644
--- a/src/ode_solver/builder.rs
+++ b/src/ode_solver/builder.rs
@@ -1,4 +1,4 @@
-use crate::{op::Op, Matrix, OdeSolverProblem, Vector};
+use crate::{op::Op, vector::DefaultDenseMatrix, Matrix, OdeSolverProblem, Vector};
 use anyhow::Result;
 
 use super::equations::{OdeSolverEquations, OdeSolverEquationsMassI};
@@ -245,6 +245,23 @@ impl OdeBuilder {
         ))
     }
 
+    /// Build an ODE problem using the default dense matrix (see [Self::build_ode]).
+    #[allow(clippy::type_complexity)]
+    pub fn build_ode_dense<V, F, G, I>(
+        self,
+        rhs: F,
+        rhs_jac: G,
+        init: I,
+    ) -> Result<OdeSolverProblem<OdeSolverEquationsMassI<V::M, F, G, I>>>
+    where
+        V: Vector + DefaultDenseMatrix,
+        F: Fn(&V, &V, V::T, &mut V),
+        G: Fn(&V, &V, V::T, &V, &mut V),
+        I: Fn(&V, V::T) -> V,
+    {
+        self.build_ode(rhs, rhs_jac, init)
+    }
+
     /// Build an ODE problem using the DiffSL language (requires the `diffsl` feature).
     /// The source code is provided as a string, please see the [DiffSL documentation](https://martinjrobins.github.io/diffsl/) for more information.
     #[cfg(feature = "diffsl")]
diff --git a/src/ode_solver/method.rs b/src/ode_solver/method.rs
index d9946cf8..d8af5603 100644
--- a/src/ode_solver/method.rs
+++ b/src/ode_solver/method.rs
@@ -3,7 +3,7 @@ use std::rc::Rc;
 
 use crate::{
     op::{filter::FilterCallable, ode_rhs::OdeRhs},
-    Matrix, NonLinearSolver, OdeEquations, OdeSolverProblem, SolverProblem, Vector, VectorIndex,
+    NonLinearSolver, OdeEquations, OdeSolverProblem, SolverProblem, Vector, VectorIndex,
 };
 
 /// Trait for ODE solver methods. This is the main user interface for the ODE solvers.
@@ -31,7 +31,7 @@ pub trait OdeSolverMethod<Eqn: OdeEquations> {
 
     /// Set the problem to solve, this performs any initialisation required by the solver. Call this before calling `step` or `solve`.
     /// The solver takes ownership of the initial state given by `state`, this is assumed to be consistent with any algebraic constraints.
-    fn set_problem(&mut self, state: OdeSolverState<Eqn::M>, problem: &OdeSolverProblem<Eqn>);
+    fn set_problem(&mut self, state: OdeSolverState<Eqn::V>, problem: &OdeSolverProblem<Eqn>);
 
     /// Step the solution forward by one step, altering the internal state of the solver.
     fn step(&mut self) -> Result<()>;
@@ -40,12 +40,12 @@ pub trait OdeSolverMethod<Eqn: OdeEquations> {
     fn interpolate(&self, t: Eqn::T) -> Result<Eqn::V>;
 
     /// Get the current state of the solver, if it exists
-    fn state(&self) -> Option<&OdeSolverState<Eqn::M>>;
+    fn state(&self) -> Option<&OdeSolverState<Eqn::V>>;
 
     /// Take the current state of the solver, if it exists, returning it to the user. This is useful if you want to use this
     /// state in another solver or problem. Note that this will unset the current problem and solver state, so you will need to call
     /// `set_problem` again before calling `step` or `solve`.
-    fn take_state(&mut self) -> Option<OdeSolverState<Eqn::M>>;
+    fn take_state(&mut self) -> Option<OdeSolverState<Eqn::V>>;
 
     /// Reinitialise the solver state and solve the problem up to time `t`
     fn solve(&mut self, problem: &OdeSolverProblem<Eqn>, t: Eqn::T) -> Result<Eqn::V> {
@@ -75,29 +75,23 @@ pub trait OdeSolverMethod<Eqn: OdeEquations> {
 
 /// State for the ODE solver, containing the current solution `y`, the current time `t`, and the current step size `h`.
 #[derive(Clone)]
-pub struct OdeSolverState<M: Matrix> {
-    pub y: M::V,
-    pub t: M::T,
-    pub h: M::T,
-    _phantom: std::marker::PhantomData<M>,
+pub struct OdeSolverState<V: Vector> {
+    pub y: V,
+    pub t: V::T,
+    pub h: V::T,
 }
 
-impl<M: Matrix> OdeSolverState<M> {
+impl<V: Vector> OdeSolverState<V> {
     /// Create a new solver state from an ODE problem. Note that this does not make the state consistent with the algebraic constraints.
     /// If you need to make the state consistent, use `new_consistent` instead.
     pub fn new<Eqn>(ode_problem: &OdeSolverProblem<Eqn>) -> Self
     where
-        Eqn: OdeEquations<M = M, T = M::T, V = M::V>,
+        Eqn: OdeEquations<T = V::T, V = V>,
     {
         let t = ode_problem.t0;
         let h = ode_problem.h0;
         let y = ode_problem.eqn.init(t);
-        Self {
-            y,
-            t,
-            h,
-            _phantom: std::marker::PhantomData,
-        }
+        Self { y, t, h }
     }
 
     /// Create a new solver state from an ODE problem, making the state consistent with the algebraic constraints.
@@ -106,7 +100,7 @@ impl<M: Matrix> OdeSolverState<M> {
         root_solver: &mut S,
     ) -> Result<Self>
     where
-        Eqn: OdeEquations<M = M, T = M::T, V = M::V>,
+        Eqn: OdeEquations<T = V::T, V = V>,
         S: NonLinearSolver<FilterCallable<OdeRhs<Eqn>>> + ?Sized,
     {
         let t = ode_problem.t0;
@@ -114,12 +108,7 @@ impl<M: Matrix> OdeSolverState<M> {
         let indices = ode_problem.eqn.algebraic_indices();
         let mut y = ode_problem.eqn.init(t);
         if indices.len() == 0 {
-            return Ok(Self {
-                y,
-                t,
-                h,
-                _phantom: std::marker::PhantomData,
-            });
+            return Ok(Self { y, t, h });
         }
         let mut y_filtered = y.filter(&indices);
         let atol = Rc::new(ode_problem.atol.as_ref().filter(&indices));
@@ -132,11 +121,6 @@ impl<M: Matrix> OdeSolverState<M> {
         let init_problem = root_solver.problem().unwrap();
         let indices = init_problem.f.indices();
         y.scatter_from(&y_filtered, indices);
-        Ok(Self {
-            y,
-            t,
-            h,
-            _phantom: std::marker::PhantomData,
-        })
+        Ok(Self { y, t, h })
     }
 }
diff --git a/src/ode_solver/sdirk.rs b/src/ode_solver/sdirk.rs
index 56f2dd66..533ec7c9 100644
--- a/src/ode_solver/sdirk.rs
+++ b/src/ode_solver/sdirk.rs
@@ -36,7 +36,7 @@ where
     tableau: Tableau<M>,
     problem: Option<OdeSolverProblem<Eqn>>,
     nonlinear_solver: NewtonNonlinearSolver<SdirkCallable<Eqn>>,
-    state: Option<OdeSolverState<Eqn::M>>,
+    state: Option<OdeSolverState<Eqn::V>>,
     diff: M,
     gamma: Eqn::T,
     is_sdirk: bool,
@@ -175,7 +175,7 @@ where
 
     fn set_problem(
         &mut self,
-        mut state: OdeSolverState<<Eqn>::M>,
+        mut state: OdeSolverState<<Eqn>::V>,
         problem: &OdeSolverProblem<Eqn>,
     ) {
         // update initial step size based on function
@@ -490,11 +490,11 @@ where
         }
     }
 
-    fn state(&self) -> Option<&OdeSolverState<<Eqn>::M>> {
+    fn state(&self) -> Option<&OdeSolverState<<Eqn>::V>> {
         self.state.as_ref()
     }
 
-    fn take_state(&mut self) -> Option<OdeSolverState<<Eqn>::M>> {
+    fn take_state(&mut self) -> Option<OdeSolverState<<Eqn>::V>> {
         Option::take(&mut self.state)
     }
 }
diff --git a/src/ode_solver/sundials.rs b/src/ode_solver/sundials.rs
index 2aad73ce..8f7f43e7 100644
--- a/src/ode_solver/sundials.rs
+++ b/src/ode_solver/sundials.rs
@@ -123,7 +123,7 @@ where
     yp: SundialsVector,
     jacobian: SundialsMatrix,
     statistics: SundialsStatistics,
-    state: Option<OdeSolverState<Eqn::M>>,
+    state: Option<OdeSolverState<Eqn::V>>,
 }
 
 impl<Eqn> SundialsIda<Eqn>
@@ -251,15 +251,15 @@ where
         self.problem.as_ref()
     }
 
-    fn state(&self) -> Option<&OdeSolverState<<Eqn>::M>> {
+    fn state(&self) -> Option<&OdeSolverState<<Eqn>::V>> {
         self.state.as_ref()
     }
 
-    fn take_state(&mut self) -> Option<OdeSolverState<<Eqn>::M>> {
+    fn take_state(&mut self) -> Option<OdeSolverState<<Eqn>::V>> {
         Option::take(&mut self.state)
     }
 
-    fn set_problem(&mut self, state: OdeSolverState<Eqn::M>, problem: &OdeSolverProblem<Eqn>) {
+    fn set_problem(&mut self, state: OdeSolverState<Eqn::V>, problem: &OdeSolverProblem<Eqn>) {
         self.state = Some(state);
         let state = self.state.as_ref().unwrap();
         self.problem = Some(problem.clone());
diff --git a/src/vector/faer_serial.rs b/src/vector/faer_serial.rs
index 6739c1b6..b5d00b23 100644
--- a/src/vector/faer_serial.rs
+++ b/src/vector/faer_serial.rs
@@ -1,11 +1,17 @@
 use std::ops::{Div, Mul, MulAssign};
 
-use faer::{unzipped, zipped, Col, ColMut, ColRef};
+use faer::{unzipped, zipped, Col, ColMut, ColRef, Mat};
 
 use crate::{scalar::Scale, IndexType, Scalar};
 
 use crate::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut};
 
+use super::DefaultDenseMatrix;
+
+impl<T: Scalar> DefaultDenseMatrix for Col<T> {
+    type M = Mat<T>;
+}
+
 macro_rules! impl_op_for_faer_struct {
     ($struct:ident, $trait_name:ident, $func_name:ident) => {
         impl<'a, T: Scalar> $trait_name<Scale<T>> for $struct<'a, T> {
diff --git a/src/vector/mod.rs b/src/vector/mod.rs
index a2a09b6e..51219703 100644
--- a/src/vector/mod.rs
+++ b/src/vector/mod.rs
@@ -1,3 +1,4 @@
+use crate::matrix::DenseMatrix;
 use crate::scalar::Scale;
 use crate::{IndexType, Scalar};
 use num_traits::Zero;
@@ -206,3 +207,7 @@ pub trait Vector:
         }
     }
 }
+
+pub trait DefaultDenseMatrix: Vector {
+    type M: DenseMatrix<V = Self, T = Self::T>;
+}
diff --git a/src/vector/nalgebra_serial.rs b/src/vector/nalgebra_serial.rs
index de5de1ca..681fd78c 100644
--- a/src/vector/nalgebra_serial.rs
+++ b/src/vector/nalgebra_serial.rs
@@ -1,10 +1,14 @@
 use std::ops::{Div, Mul, MulAssign};
 
-use nalgebra::{DVector, DVectorView, DVectorViewMut};
+use nalgebra::{DMatrix, DVector, DVectorView, DVectorViewMut};
 
 use crate::{IndexType, Scalar, Scale};
 
-use super::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut};
+use super::{DefaultDenseMatrix, Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut};
+
+impl<T: Scalar> DefaultDenseMatrix for DVector<T> {
+    type M = DMatrix<T>;
+}
 
 macro_rules! impl_op_for_dvector_struct {
     ($struct:ident, $trait_name:ident, $func_name:ident) => {