diff --git a/src/jacobian/mod.rs b/src/jacobian/mod.rs
index f8a7a728..77ac42f8 100644
--- a/src/jacobian/mod.rs
+++ b/src/jacobian/mod.rs
@@ -46,7 +46,7 @@ pub fn find_non_zeros_linear<F: LinearOp + ?Sized>(op: &F, t: F::T) -> Vec<(usiz
     let mut triplets = Vec::with_capacity(op.nstates());
     for j in 0..op.nstates() {
         v[j] = F::T::NAN;
-        op.call_inplace(&v, t, &mut col);
+        op.call_inplace(v.view(), t, col.view_mut());
         for i in 0..op.nout() {
             if col[i].is_nan() {
                 triplets.push((i, j));
@@ -144,7 +144,7 @@ impl<M: Matrix> JacobianColoring<M> {
             let dst_indices = &self.dst_indices_per_color[c];
             let src_indices = &self.src_indices_per_color[c];
             v.assign_at_indices(input, F::T::one());
-            op.call_inplace(&v, t, &mut col);
+            op.call_inplace(v.view(), t, col.view_mut());
             y.set_data_with_indices(dst_indices, src_indices, &col);
             v.assign_at_indices(input, F::T::zero());
         }
@@ -165,7 +165,7 @@ mod tests {
         jacobian::{coloring::nonzeros2graph, greedy_coloring::color_graph_greedy},
         op::closure::Closure,
     };
-    use nalgebra::{DMatrix, DVector};
+    use nalgebra::{DMatrix, DVector, DVectorView, DVectorViewMut};
     use std::ops::MulAssign;
 
     fn helper_triplets2op_nonlinear(
@@ -175,17 +175,17 @@ mod tests {
     ) -> impl NonLinearOp<M = DMatrix<f64>, V = DVector<f64>, T = f64> + '_ {
         let nstates = ncols;
         let nout = nrows;
-        let f = move |x: &DVector<f64>, y: &mut DVector<f64>| {
+        let f = move |x: DVectorView<'_, f64>, mut y: DVectorViewMut<'_, f64>| {
             for (i, j, v) in triplets {
                 y[*i] += x[*j] * v;
             }
         };
         let mut ret = Closure::new(
-            move |x: &DVector<f64>, _p: &DVector<f64>, _t, y: &mut DVector<f64>| {
+            move |x: DVectorView<'_, f64>, _p: &DVector<f64>, _t, mut y: DVectorViewMut<'_, f64>| {
                 y.fill(0.0);
                 f(x, y);
             },
-            move |_x: &DVector<f64>, _p: &DVector<f64>, _t, v, y: &mut DVector<f64>| {
+            move |_x: DVectorView<'_, f64>, _p: &DVector<f64>, _t, v, mut y: DVectorViewMut<'_, f64>| {
                 y.fill(0.0);
                 f(v, y);
             },
@@ -302,7 +302,7 @@ mod tests {
             coloring.matrix_inplace(&op, t0, &mut jac);
             let mut gemv1 = V::zeros(n);
             let v = V::from_element(3, 1.0);
-            op.gemv_inplace(&v, t0, 0.0, &mut gemv1);
+            op.gemv_inplace(Vector::view(&v), t0, 0.0, Vector::view_mut(&mut gemv1));
             let mut gemv2 = V::zeros(n);
             jac.gemv(1.0, &v, 0.0, &mut gemv2);
             gemv1.assert_eq_st(&gemv2, 1e-10);
diff --git a/src/matrix/dense_faer_serial.rs b/src/matrix/dense_faer_serial.rs
index 7e8d0e89..9ca776ee 100644
--- a/src/matrix/dense_faer_serial.rs
+++ b/src/matrix/dense_faer_serial.rs
@@ -2,7 +2,7 @@ use std::ops::{Mul, MulAssign};
 
 use super::default_solver::DefaultSolver;
 use super::{Dense, DenseMatrix, Matrix, MatrixCommon, MatrixSparsity, MatrixView, MatrixViewMut};
-use crate::op::NonLinearOp;
+use crate::op::{NonLinearOp, VView, VViewMut};
 use crate::scalar::{IndexType, Scalar, Scale};
 use crate::vector::Vector;
 use crate::FaerLU;
@@ -155,7 +155,7 @@ impl<T: Scalar> Matrix for Mat<T> {
         }
         Ok(m)
     }
-    fn gemv(&self, alpha: Self::T, x: &Self::V, beta: Self::T, y: &mut Self::V) {
+    fn gemv(&self, alpha: Self::T, x: VView<'_, Self>, beta: Self::T, mut y: VViewMut<'_, Self>) {
         *y = faer::scale(alpha) * self * x + faer::scale(beta) * &*y;
     }
     fn zeros(nrows: IndexType, ncols: IndexType) -> Self {
diff --git a/src/matrix/dense_nalgebra_serial.rs b/src/matrix/dense_nalgebra_serial.rs
index 61e0da3c..856d435c 100644
--- a/src/matrix/dense_nalgebra_serial.rs
+++ b/src/matrix/dense_nalgebra_serial.rs
@@ -3,7 +3,7 @@ use std::ops::{AddAssign, Mul, MulAssign};
 use anyhow::Result;
 use nalgebra::{DMatrix, DMatrixView, DMatrixViewMut, DVector, DVectorView, DVectorViewMut};
 
-use crate::op::NonLinearOp;
+use crate::op::{NonLinearOp, VView, VViewMut};
 use crate::{scalar::Scale, IndexType, Scalar};
 
 use crate::{DenseMatrix, Matrix, MatrixCommon, MatrixView, MatrixViewMut, NalgebraLU};
@@ -127,8 +127,8 @@ impl<T: Scalar> Matrix for DMatrix<T> {
         self.diagonal()
     }
 
-    fn gemv(&self, alpha: Self::T, x: &Self::V, beta: Self::T, y: &mut Self::V) {
-        y.gemv(alpha, self, x, beta);
+    fn gemv(&self, alpha: Self::T, x: VView<'_, Self>, beta: Self::T, mut y: VViewMut<'_, Self>) {
+        y.gemv(alpha, self, &x, beta);
     }
     fn copy_from(&mut self, other: &Self) {
         self.copy_from(other);
diff --git a/src/matrix/mod.rs b/src/matrix/mod.rs
index aa8ef137..ecdd6ae2 100644
--- a/src/matrix/mod.rs
+++ b/src/matrix/mod.rs
@@ -1,6 +1,7 @@
 use std::fmt::Debug;
 use std::ops::{Add, AddAssign, Index, IndexMut, Mul, MulAssign, Sub, SubAssign};
 
+use crate::op::{VView, VViewMut};
 use crate::scalar::Scale;
 use crate::{IndexType, Scalar, Vector};
 use anyhow::Result;
@@ -210,7 +211,7 @@ pub trait Matrix:
     fn diagonal(&self) -> Self::V;
 
     /// Perform a matrix-vector multiplication `y = self * x + beta * y`.
-    fn gemv(&self, alpha: Self::T, x: &Self::V, beta: Self::T, y: &mut Self::V);
+    fn gemv(&self, alpha: Self::T, x: VView<'_, Self>, beta: Self::T, y: VViewMut<'_, Self>);
 
     /// Copy the contents of `other` into `self`
     fn copy_from(&mut self, other: &Self);
diff --git a/src/matrix/sparse_serial.rs b/src/matrix/sparse_serial.rs
index 13d5dcf2..93968831 100644
--- a/src/matrix/sparse_serial.rs
+++ b/src/matrix/sparse_serial.rs
@@ -4,7 +4,7 @@ use anyhow::Result;
 use nalgebra::DVector;
 use nalgebra_sparse::{pattern::SparsityPattern, CooMatrix, CscMatrix};
 
-use crate::{scalar::Scale, IndexType, Scalar};
+use crate::{op::{VView, VViewMut}, scalar::Scale, IndexType, Scalar};
 
 use super::{Matrix, MatrixCommon, MatrixSparsity};
 
@@ -180,7 +180,7 @@ impl<T: Scalar> Matrix for CscMatrix<T> {
     fn copy_from(&mut self, other: &Self) {
         self.clone_from(other);
     }
-    fn gemv(&self, alpha: Self::T, x: &Self::V, beta: Self::T, y: &mut Self::V) {
+    fn gemv(&self, alpha: Self::T, x: VView<'_, Self>, beta: Self::T, mut y: VViewMut<'_, Self>) {
         let mut tmp = self * x;
         tmp *= alpha;
         y.axpy(alpha, &tmp, beta);
diff --git a/src/matrix/sundials.rs b/src/matrix/sundials.rs
index 0c80c788..91093aae 100644
--- a/src/matrix/sundials.rs
+++ b/src/matrix/sundials.rs
@@ -12,7 +12,7 @@ use sundials_sys::{
 
 use crate::{
     ode_solver::sundials::sundials_check,
-    op::NonLinearOp,
+    op::{NonLinearOp, VView, VViewMut},
     scalar::scale,
     vector::sundials::{get_suncontext, SundialsVector},
     IndexType, Scale, SundialsLinearSolver, Vector,
@@ -302,7 +302,7 @@ impl Matrix for SundialsMatrix {
     }
 
     /// Perform a matrix-vector multiplication `y = alpha * self * x + beta * y`.
-    fn gemv(&self, alpha: Self::T, x: &Self::V, beta: Self::T, y: &mut Self::V) {
+    fn gemv(&self, alpha: Self::T, x: VView<'_, Self>, beta: Self::T, mut y: VViewMut<'_, Self>) {
         let a = self.sundials_matrix();
         let tmp = SundialsVector::new_serial(self.nrows());
         sundials_check(unsafe { SUNMatMatvecSetup(a) }).unwrap();
diff --git a/src/nonlinear_solver/newton.rs b/src/nonlinear_solver/newton.rs
index 86312f35..783fe683 100644
--- a/src/nonlinear_solver/newton.rs
+++ b/src/nonlinear_solver/newton.rs
@@ -74,7 +74,7 @@ impl<C: NonLinearOp, Ls: LinearSolver<C>> NonLinearSolver<C> for NewtonNonlinear
         self.niter = 0;
         loop {
             self.niter += 1;
-            problem.f.call_inplace(xn, t, &mut tmp);
+            problem.f.call_inplace(xn.view(), t, tmp.view_mut());
             //tmp = f_at_n
 
             self.linear_solver.solve_in_place(&mut tmp)?;
diff --git a/src/nonlinear_solver/root.rs b/src/nonlinear_solver/root.rs
index 88f74448..21140f3b 100644
--- a/src/nonlinear_solver/root.rs
+++ b/src/nonlinear_solver/root.rs
@@ -27,7 +27,7 @@ impl<V: Vector> RootFinder<V> {
     /// Set the lower boundary of the root search.
     /// This function should be called first after [Self::new]
     pub fn init(&self, root_fn: &impl NonLinearOp<V = V, T = V::T>, y: &V, t: V::T) {
-        root_fn.call_inplace(y, t, &mut self.g0.borrow_mut());
+        root_fn.call_inplace(y.view(), t, self.g0.borrow_mut().view_mut());
         self.t0.replace(t);
     }
 
@@ -48,7 +48,7 @@ impl<V: Vector> RootFinder<V> {
         let g1 = &mut *self.g1.borrow_mut();
         let g0 = &mut *self.g0.borrow_mut();
         let gmid = &mut *self.gmid.borrow_mut();
-        root_fn.call_inplace(y, t, g1);
+        root_fn.call_inplace(y.view(), t, g1.view_mut());
 
         let sign_change_fn = |mut acc: (bool, V::T, i32), g0: V::T, g1: V::T, i: IndexType| {
             if g1 == V::T::zero() {
@@ -115,7 +115,7 @@ impl<V: Vector> RootFinder<V> {
             }
 
             let ymid = interpolate(t_mid).unwrap();
-            root_fn.call_inplace(&ymid, t_mid, gmid);
+            root_fn.call_inplace(ymid.view(), t_mid, gmid.view_mut());
 
             let (rootfnd, _gfracmax, imax_i32) =
                 (*g0).binary_fold(gmid, (false, V::T::zero(), -1), sign_change_fn);
@@ -128,7 +128,7 @@ impl<V: Vector> RootFinder<V> {
                 std::mem::swap(g1, gmid);
             } else if rootfnd {
                 // we are returning so make sure g0 is set for next iteration
-                root_fn.call_inplace(y, t, g0);
+                root_fn.call_inplace(y.view(), t, g0.view_mut());
 
                 // No sign change in (tlo,tmid), but g = 0 at tmid; return root tmid.
                 return Some(t_mid);
@@ -151,7 +151,7 @@ impl<V: Vector> RootFinder<V> {
             i += 1;
         }
         // we are returning so make sure g0 is set for next iteration
-        root_fn.call_inplace(y, t, g0);
+        root_fn.call_inplace(y.view(), t, g0.view_mut());
         Some(t1)
     }
 }
diff --git a/src/ode_solver/bdf/mod.rs b/src/ode_solver/bdf/mod.rs
index 0a53ebb5..0fdb7f8f 100644
--- a/src/ode_solver/bdf/mod.rs
+++ b/src/ode_solver/bdf/mod.rs
@@ -358,11 +358,11 @@ where
         scale_factor *= scale(problem.rtol);
         scale_factor += problem.atol.as_ref();
 
-        let f0 = problem.eqn.rhs().call(&state.y, state.t);
+        let f0 = problem.eqn.rhs().call(state.y.view(), state.t);
         let hf0 = &f0 * scale(state.h);
         let y1 = &state.y + &hf0;
         let t1 = state.t + state.h;
-        let f1 = problem.eqn.rhs().call(&y1, t1);
+        let f1 = problem.eqn.rhs().call(y1.view(), t1);
 
         // store f1 in diff[1] for use in step size control
         self.diff.column_mut(1).copy_from(&hf0);
diff --git a/src/ode_solver/diffsl.rs b/src/ode_solver/diffsl.rs
index 34dbdca1..def11951 100644
--- a/src/ode_solver/diffsl.rs
+++ b/src/ode_solver/diffsl.rs
@@ -4,9 +4,7 @@ use anyhow::Result;
 use diffsl::execution::Compiler;
 
 use crate::{
-    jacobian::{find_non_zeros_linear, find_non_zeros_nonlinear, JacobianColoring},
-    op::{LinearOp, NonLinearOp, Op},
-    OdeEquations,
+    jacobian::{find_non_zeros_linear, find_non_zeros_nonlinear, JacobianColoring}, op::{LinearOp, NonLinearOp, Op}, vector::Vector, OdeEquations
 };
 
 pub type T = f64;
@@ -170,7 +168,7 @@ impl Op for DiffSlRoot<'_> {
 }
 
 impl NonLinearOp for DiffSlRoot<'_> {
-    fn call_inplace(&self, x: &Self::V, t: Self::T, y: &mut Self::V) {
+    fn call_inplace(&self, x: <Self::V as Vector>::View<'_>, t: Self::T, y: <Self::V as Vector>::ViewMut<'_>) {
         self.context.compiler.calc_stop(
             t,
             x.as_slice(),
@@ -185,7 +183,7 @@ impl NonLinearOp for DiffSlRoot<'_> {
 }
 
 impl NonLinearOp for DiffSlRhs<'_> {
-    fn call_inplace(&self, x: &Self::V, t: Self::T, y: &mut Self::V) {
+    fn call_inplace(&self, x: <Self::V as Vector>::View<'_>, t: Self::T, y: <Self::V as Vector>::ViewMut<'_>) {
         self.context.compiler.rhs(
             t,
             x.as_slice(),
@@ -217,7 +215,7 @@ impl NonLinearOp for DiffSlRhs<'_> {
 }
 
 impl LinearOp for DiffSlMass<'_> {
-    fn gemv_inplace(&self, x: &Self::V, t: Self::T, beta: Self::T, y: &mut Self::V) {
+    fn gemv_inplace(&self, x: <Self::V as Vector>::View<'_>, t: Self::T, beta: Self::T, y: <Self::V as Vector>::ViewMut<'_>) {
         let mut tmp = self.context.tmp.borrow_mut();
         self.context.compiler.mass(
             t,
@@ -339,7 +337,7 @@ mod tests {
         let init = eqn.init(0.0);
         let init_expect = DVector::from_vec(vec![y0, 0.0]);
         init.assert_eq_st(&init_expect, 1e-10);
-        let rhs = eqn.rhs().call(&init, 0.0);
+        let rhs = eqn.rhs().call(Vector::view(&init), 0.0);
         let rhs_expect = DVector::from_vec(vec![r * y0 * (1.0 - y0 / k), 2.0 * y0]);
         rhs.assert_eq_st(&rhs_expect, 1e-10);
         let v = DVector::from_vec(vec![1.0, 1.0]);
@@ -348,7 +346,7 @@ mod tests {
         rhs_jac.assert_eq_st(&rhs_jac_expect, 1e-10);
         let mut mass_y = DVector::from_vec(vec![0.0, 0.0]);
         let v = DVector::from_vec(vec![1.0, 1.0]);
-        eqn.mass().call_inplace(&v, 0.0, &mut mass_y);
+        eqn.mass().call_inplace(Vector::view(&v), 0.0, Vector::view_mut(&mut mass_y));
         let mass_y_expect = DVector::from_vec(vec![1.0, 0.0]);
         mass_y.assert_eq_st(&mass_y_expect, 1e-10);
 
diff --git a/src/ode_solver/equations.rs b/src/ode_solver/equations.rs
index 42e39dfc..27baa3c5 100644
--- a/src/ode_solver/equations.rs
+++ b/src/ode_solver/equations.rs
@@ -181,7 +181,7 @@ mod tests {
     fn ode_equation_test() {
         let (problem, _soln) = exponential_decay_problem::<Mcpu>(false);
         let y = DVector::from_vec(vec![1.0, 1.0]);
-        let rhs_y = problem.eqn.rhs().call(&y, 0.0);
+        let rhs_y = problem.eqn.rhs().call(Vector::view(&y), 0.0);
         let expect_rhs_y = DVector::from_vec(vec![-0.1, -0.1]);
         rhs_y.assert_eq_st(&expect_rhs_y, 1e-10);
         let jac_rhs_y = problem.eqn.rhs().jac_mul(&y, 0.0, &y);
@@ -203,7 +203,7 @@ mod tests {
     fn ode_with_mass_test() {
         let (problem, _soln) = exponential_decay_with_algebraic_problem::<Mcpu>(false);
         let y = DVector::from_vec(vec![1.0, 1.0, 1.0]);
-        let rhs_y = problem.eqn.rhs().call(&y, 0.0);
+        let rhs_y = problem.eqn.rhs().call(Vector::view(&y), 0.0);
         let expect_rhs_y = DVector::from_vec(vec![-0.1, -0.1, 0.0]);
         rhs_y.assert_eq_st(&expect_rhs_y, 1e-10);
         let jac_rhs_y = problem.eqn.rhs().jac_mul(&y, 0.0, &y);
diff --git a/src/ode_solver/mod.rs b/src/ode_solver/mod.rs
index 178d5974..848c53c1 100644
--- a/src/ode_solver/mod.rs
+++ b/src/ode_solver/mod.rs
@@ -533,7 +533,7 @@ mod tests {
     }
 
     impl<M: Matrix> NonLinearOp for TestEqnRhs<M> {
-        fn call_inplace(&self, _x: &Self::V, _t: Self::T, y: &mut Self::V) {
+        fn call_inplace(&self, _x: <Self::V as Vector>::View<'_>, _t: Self::T, y: <Self::V as Vector>::ViewMut<'_>) {
             y[0] = M::T::zero();
         }
 
diff --git a/src/ode_solver/sdirk.rs b/src/ode_solver/sdirk.rs
index 4cdd34d4..3c1e3478 100644
--- a/src/ode_solver/sdirk.rs
+++ b/src/ode_solver/sdirk.rs
@@ -218,7 +218,7 @@ where
         // compute first step based on alg in Hairer, Norsett, Wanner
         // Solving Ordinary Differential Equations I, Nonstiff Problems
         // Section II.4.2
-        let f0 = problem.eqn.rhs().call(&state.y, state.t);
+        let f0 = problem.eqn.rhs().call(state.y.view(), state.t);
         let hf0 = &f0 * scale(state.h);
 
         let mut tmp = f0.clone();
@@ -237,7 +237,7 @@ where
 
         let y1 = &state.y + hf0;
         let t1 = state.t + h0;
-        let f1 = problem.eqn.rhs().call(&y1, t1);
+        let f1 = problem.eqn.rhs().call(y1.view(), t1);
 
         let mut df = f1 - &f0;
         df *= scale(Eqn::T::one() / h0);
diff --git a/src/ode_solver/sundials.rs b/src/ode_solver/sundials.rs
index 0ad09284..1ceb90f1 100644
--- a/src/ode_solver/sundials.rs
+++ b/src/ode_solver/sundials.rs
@@ -153,9 +153,9 @@ where
         let mut rr = SundialsVector::new_not_owned(rr);
         // F(t, y, y') =  M y' - f(t, y)
         // rr = f(t, y)
-        data.eqn.rhs().call_inplace(&y, t, &mut rr);
+        data.eqn.rhs().call_inplace(y.view(), t, rr.view_mut());
         // rr = M y' - rr
-        data.eqn.mass().gemv_inplace(&yp, t, -1.0, &mut rr);
+        data.eqn.mass().gemv_inplace(yp.view(), t, -1.0, rr.view_mut());
         0
     }
 
diff --git a/src/op/bdf.rs b/src/op/bdf.rs
index 698d2f5f..cb6626ad 100644
--- a/src/op/bdf.rs
+++ b/src/op/bdf.rs
@@ -120,25 +120,25 @@ where
     for<'b> &'b Eqn::M: MatrixRef<Eqn::M>,
 {
     // F(y) = M (y - y0 + psi) - c * f(y) = 0
-    fn call_inplace(&self, x: &Eqn::V, t: Eqn::T, y: &mut Eqn::V) {
+    fn call_inplace(&self, x: <Eqn::V as Vector>::View<'_>, t: Eqn::T, y: <Eqn::V as Vector>::ViewMut<'_>) {
         let psi_neg_y0_ref = self.psi_neg_y0.borrow();
         let psi_neg_y0 = psi_neg_y0_ref.deref();
 
         self.eqn.rhs().call_inplace(x, t, y);
 
         let mut tmp = self.tmp.borrow_mut();
-        tmp.copy_from(x);
+        tmp.copy_from_view(&x);
         tmp.add_assign(psi_neg_y0);
         let c = *self.c.borrow().deref();
         // y = M tmp - c * y
-        self.eqn.mass().gemv_inplace(&tmp, t, -c, y);
+        self.eqn.mass().gemv_inplace(tmp.view(), t, -c, y);
     }
     // (M - c * f'(y)) v
     fn jac_mul_inplace(&self, x: &Eqn::V, t: Eqn::T, v: &Eqn::V, y: &mut Eqn::V) {
         self.eqn.rhs().jac_mul_inplace(x, t, v, y);
         let c = *self.c.borrow().deref();
         // y = Mv - c y
-        self.eqn.mass().gemv_inplace(v, t, -c, y);
+        self.eqn.mass().gemv_inplace(v.view(), t, -c, y.view_mut());
     }
 
     fn jacobian_inplace(&self, x: &Self::V, t: Self::T, y: &mut Self::M) {
@@ -200,7 +200,7 @@ mod tests {
         //        |-0.1|
         //  i.e. F(y) = |1 0| |2.1| - 0.1 * |-0.1| =  |2.11|
         //              |0 1| |2.2|         |-0.1|    |2.21|
-        bdf_callable.call_inplace(&y, t, &mut y_out);
+        bdf_callable.call_inplace(Vector::view(&y), t, Vector::view_mut(&mut y_out));
         let y_out_expect = Vcpu::from_vec(vec![2.11, 2.21]);
         y_out.assert_eq_st(&y_out_expect, 1e-10);
 
diff --git a/src/op/closure.rs b/src/op/closure.rs
index 039e199c..037da916 100644
--- a/src/op/closure.rs
+++ b/src/op/closure.rs
@@ -6,13 +6,13 @@ use crate::{
     Matrix, Vector,
 };
 
-use super::{NonLinearOp, Op, OpStatistics};
+use super::{NonLinearOp, Op, OpStatistics, VView, VViewMut};
 
 pub struct Closure<M, F, G>
 where
     M: Matrix,
-    F: Fn(&M::V, &M::V, M::T, &mut M::V),
-    G: Fn(&M::V, &M::V, M::T, &M::V, &mut M::V),
+    F: Fn(VView<'_, M>, &M::V, M::T, VViewMut<'_, M>),
+    G: Fn(VView<'_, M>, &M::V, M::T, VView<'_, M>, VViewMut<'_, M>),
 {
     func: F,
     jacobian_action: G,
@@ -28,8 +28,8 @@ where
 impl<M, F, G> Closure<M, F, G>
 where
     M: Matrix,
-    F: Fn(&M::V, &M::V, M::T, &mut M::V),
-    G: Fn(&M::V, &M::V, M::T, &M::V, &mut M::V),
+    F: Fn(VView<'_, M>, &M::V, M::T, VViewMut<'_, M>),
+    G: Fn(VView<'_, M>, &M::V, M::T, VView<'_, M>, VViewMut<'_, M>),
 {
     pub fn new(func: F, jacobian_action: G, nstates: usize, nout: usize, p: Rc<M::V>) -> Self {
         let nparams = p.len();
@@ -59,8 +59,8 @@ where
 impl<M, F, G> Op for Closure<M, F, G>
 where
     M: Matrix,
-    F: Fn(&M::V, &M::V, M::T, &mut M::V),
-    G: Fn(&M::V, &M::V, M::T, &M::V, &mut M::V),
+    F: Fn(<M::V as Vector>::View<'_>, &M::V, M::T, <M::V as Vector>::ViewMut<'_>),
+    G: Fn(<M::V as Vector>::View<'_>, &M::V, M::T, <M::V as Vector>::View<'_>, <M::V as Vector>::ViewMut<'_>),
 {
     type V = M::V;
     type T = M::T;
@@ -89,10 +89,10 @@ where
 impl<M, F, G> NonLinearOp for Closure<M, F, G>
 where
     M: Matrix,
-    F: Fn(&M::V, &M::V, M::T, &mut M::V),
-    G: Fn(&M::V, &M::V, M::T, &M::V, &mut M::V),
+    F: Fn(VView<'_, M>, &M::V, M::T, VViewMut<'_, M>),
+    G: Fn(VView<'_, M>, &M::V, M::T, VView<'_, M>, VViewMut<'_, M>),
 {
-    fn call_inplace(&self, x: &M::V, t: M::T, y: &mut M::V) {
+    fn call_inplace(&self, x: <M::V as Vector>::View<'_>, t: M::T, y: <M::V as Vector>::ViewMut<'_>) {
         self.statistics.borrow_mut().increment_call();
         (self.func)(x, self.p.as_ref(), t, y)
     }
diff --git a/src/op/mod.rs b/src/op/mod.rs
index ef260a37..0d602e76 100644
--- a/src/op/mod.rs
+++ b/src/op/mod.rs
@@ -1,6 +1,6 @@
 use std::rc::Rc;
 
-use crate::{Matrix, Scalar, Vector};
+use crate::{matrix::MatrixCommon, Matrix, Scalar, Vector};
 
 use num_traits::{One, Zero};
 use serde::Serialize;
@@ -16,6 +16,9 @@ pub mod matrix;
 pub mod sdirk;
 pub mod unit;
 
+pub type VView<'a, M> = <<M as MatrixCommon>::V as Vector>::View<'a>;
+pub type VViewMut<'a, M> = <<M as MatrixCommon>::V as Vector>::ViewMut<'a>;
+
 /// Op is a trait for operators that, given a paramter vector `p`, operates on an input vector `x` to produce an output vector `y`.
 /// It defines the number of states (i.e. length of `x`), the number of outputs (i.e. length of `y`), and number of parameters (i.e. length of `p`) of the operator.
 /// It also defines the type of the scalar, vector, and matrices used in the operator.
@@ -86,15 +89,15 @@ impl OpStatistics {
 // jac_mul_inplace method, which computes the product of the Jacobian with a given vector.
 pub trait NonLinearOp: Op {
     /// Compute the operator at a given state and time.
-    fn call_inplace(&self, x: &Self::V, t: Self::T, y: &mut Self::V);
+    fn call_inplace(&self, x: <Self::V as Vector>::View<'_>, t: Self::T, y: <Self::V as Vector>::ViewMut<'_>);
 
     /// Compute the product of the Jacobian with a given vector.
     fn jac_mul_inplace(&self, x: &Self::V, t: Self::T, v: &Self::V, y: &mut Self::V);
 
     /// Compute the operator at a given state and time, and return the result.
-    fn call(&self, x: &Self::V, t: Self::T) -> Self::V {
+    fn call(&self, x: <Self::V as Vector>::View<'_>, t: Self::T) -> Self::V {
         let mut y = Self::V::zeros(self.nout());
-        self.call_inplace(x, t, &mut y);
+        self.call_inplace(x, t, y.view_mut());
         y
     }
 
@@ -136,13 +139,13 @@ pub trait NonLinearOp: Op {
 /// calling the operator via a GEMV operation (i.e. `y = t * A * x + beta * y`), and for computing the matrix representation of the operator.
 pub trait LinearOp: Op {
     /// Compute the operator at a given state and time
-    fn call_inplace(&self, x: &Self::V, t: Self::T, y: &mut Self::V) {
+    fn call_inplace(&self, x: <Self::V as Vector>::View<'_>, t: Self::T, y: <Self::V as Vector>::ViewMut<'_>) {
         let beta = Self::T::zero();
         self.gemv_inplace(x, t, beta, y);
     }
 
     /// Computer the operator via a GEMV operation (i.e. `y = t * A * x + beta * y`)
-    fn gemv_inplace(&self, x: &Self::V, t: Self::T, beta: Self::T, y: &mut Self::V);
+    fn gemv_inplace(&self, x: <Self::V as Vector>::View<'_>, t: Self::T, beta: Self::T, y: <Self::V as Vector>::ViewMut<'_>);
 
     /// Compute the matrix representation of the operator and return it.
     fn matrix(&self, t: Self::T) -> Self::M {
@@ -162,7 +165,7 @@ pub trait LinearOp: Op {
         let mut col = Self::V::zeros(self.nout());
         for j in 0..self.nstates() {
             v[j] = Self::T::one();
-            self.call_inplace(&v, t, &mut col);
+            self.call_inplace(v.view(), t, col.view_mut());
             y.set_column(j, &col);
             v[j] = Self::T::zero();
         }
diff --git a/src/vector/faer_serial.rs b/src/vector/faer_serial.rs
index bafe7f59..0011ddcf 100644
--- a/src/vector/faer_serial.rs
+++ b/src/vector/faer_serial.rs
@@ -89,10 +89,10 @@ impl<T: Scalar> Vector for Col<T> {
     fn abs(&self) -> Self {
         zipped!(self).map(|unzipped!(xi)| xi.faer_abs())
     }
-    fn as_view(&self) -> Self::View<'_> {
+    fn view(&self) -> Self::View<'_> {
         self.as_ref()
     }
-    fn as_view_mut(&mut self) -> Self::ViewMut<'_> {
+    fn view_mut(&mut self) -> Self::ViewMut<'_> {
         self.as_mut()
     }
     fn copy_from(&mut self, other: &Self) {
@@ -123,10 +123,10 @@ impl<T: Scalar> Vector for Col<T> {
         zipped!(self).map(|unzipped!(xi)| xi.exp())
     }
     fn component_mul_assign(&mut self, other: &Self) {
-        zipped!(self.as_mut(), other.as_view()).for_each(|unzipped!(mut s, o)| *s *= *o);
+        zipped!(self.as_mut(), other.view()).for_each(|unzipped!(mut s, o)| *s *= *o);
     }
     fn component_div_assign(&mut self, other: &Self) {
-        zipped!(self.as_mut(), other.as_view()).for_each(|unzipped!(mut s, o)| *s /= *o);
+        zipped!(self.as_mut(), other.view()).for_each(|unzipped!(mut s, o)| *s /= *o);
     }
     fn filter_indices<F: Fn(Self::T) -> bool>(&self, f: F) -> Self::Index {
         let mut indices = vec![];
diff --git a/src/vector/mod.rs b/src/vector/mod.rs
index 145919a6..0b65349c 100644
--- a/src/vector/mod.rs
+++ b/src/vector/mod.rs
@@ -3,7 +3,7 @@ use crate::scalar::Scale;
 use crate::{IndexType, Scalar};
 use num_traits::Zero;
 use std::fmt::Debug;
-use std::ops::{Add, AddAssign, Div, Index, IndexMut, Mul, MulAssign, Sub, SubAssign};
+use std::ops::{Add, AddAssign, Deref, Div, Index, IndexMut, Mul, MulAssign, Sub, SubAssign};
 
 #[cfg(feature = "faer")]
 mod faer_serial;
@@ -100,6 +100,7 @@ pub trait VectorView<'a>:
     fn into_owned(self) -> Self::Owned;
 }
 
+
 pub trait Vector:
     VectorOpsByValue<Self>
     + for<'b> VectorOpsByValue<&'b Self>
@@ -136,8 +137,8 @@ pub trait Vector:
     fn zeros(nstates: usize) -> Self {
         Self::from_element(nstates, Self::T::zero())
     }
-    fn as_view(&self) -> Self::View<'_>;
-    fn as_view_mut(&mut self) -> Self::ViewMut<'_>;
+    fn view(&self) -> Self::View<'_>;
+    fn view_mut(&mut self) -> Self::ViewMut<'_>;
     fn copy_from(&mut self, other: &Self);
     fn copy_from_view(&mut self, other: &Self::View<'_>);
     fn from_vec(vec: Vec<Self::T>) -> Self;
diff --git a/src/vector/nalgebra_serial.rs b/src/vector/nalgebra_serial.rs
index 0b917ae3..d58bc4df 100644
--- a/src/vector/nalgebra_serial.rs
+++ b/src/vector/nalgebra_serial.rs
@@ -124,10 +124,10 @@ impl<T: Scalar> Vector for DVector<T> {
     fn abs(&self) -> Self {
         self.abs()
     }
-    fn as_view(&self) -> Self::View<'_> {
+    fn view(&self) -> Self::View<'_> {
         self.as_view()
     }
-    fn as_view_mut(&mut self) -> Self::ViewMut<'_> {
+    fn view_mut(&mut self) -> Self::ViewMut<'_> {
         self.as_view_mut()
     }
     fn copy_from(&mut self, other: &Self) {
diff --git a/src/vector/sundials.rs b/src/vector/sundials.rs
index ae8203cd..72fde243 100644
--- a/src/vector/sundials.rs
+++ b/src/vector/sundials.rs
@@ -84,10 +84,10 @@ impl<'a> SundialsVectorViewMut<'a> {
 pub struct SundialsVectorView<'a>(&'a SundialsVector);
 
 impl<'a> SundialsVectorView<'a> {
-    fn sundials_vector(&self) -> N_Vector {
+    pub fn sundials_vector(&self) -> N_Vector {
         self.0.sundials_vector()
     }
-    fn len(&self) -> IndexType {
+    pub fn len(&self) -> IndexType {
         unsafe { N_VGetLength_Serial(self.sundials_vector()) as IndexType }
     }
 }
@@ -458,10 +458,10 @@ impl Vector for SundialsVector {
     fn add_scalar_mut(&mut self, scalar: Self::T) {
         unsafe { N_VAddConst(self.sundials_vector(), scalar, self.sundials_vector()) }
     }
-    fn as_view(&self) -> Self::View<'_> {
+    fn view(&self) -> Self::View<'_> {
         SundialsVectorView(self)
     }
-    fn as_view_mut(&mut self) -> Self::ViewMut<'_> {
+    fn view_mut(&mut self) -> Self::ViewMut<'_> {
         SundialsVectorViewMut(self)
     }
     fn axpy(&mut self, alpha: Self::T, x: &Self, beta: Self::T) {