diff --git a/MathLib/ODE/FunctionHandles.h b/MathLib/ODE/FunctionHandles.h
index c268f21b68f043a66b3245cbf50fd16e1aa83fb8..03023aeb085d9f7eff5d66ee961b4a192be29b66 100644
--- a/MathLib/ODE/FunctionHandles.h
+++ b/MathLib/ODE/FunctionHandles.h
@@ -44,7 +44,11 @@ struct FunctionHandlesImpl : FunctionHandles
 	bool call(const double t, const double* const y,
 	          double* const ydot) override
 	{
-		if (f) return f(t, MappedConstVector<N>{y}, MappedVector<N>{ydot});
+		if (f)
+		{
+			MappedVector<N> ydot_mapped{ydot};
+			return f(t, MappedConstVector<N>{y}, ydot_mapped);
+		}
 		return false;
 	}
 
@@ -52,10 +56,13 @@ struct FunctionHandlesImpl : FunctionHandles
 	                  double* const jac) override
 	{
 		if (df)
+		{
+			MappedMatrix<N, N> jac_mapped{jac};
 			return df(t,
 			          MappedConstVector<N>{y},
 			          MappedConstVector<N>{ydot},
-			          MappedMatrix<N, N>{jac});
+			          jac_mapped);
+		}
 		return false;
 	}
 
diff --git a/MathLib/ODE/OdeSolverTypes.h b/MathLib/ODE/OdeSolverTypes.h
index 5d8affc7876e402e8eba01aedb6901521ab865db..b7d831636443138aa85f0cc2203f180ba06ff2df 100644
--- a/MathLib/ODE/OdeSolverTypes.h
+++ b/MathLib/ODE/OdeSolverTypes.h
@@ -30,13 +30,13 @@ using MappedConstVector = MappedConstMatrix<N, 1>;
 
 template <unsigned N>
 using Function = std::function<bool(
-    const double t, MappedConstVector<N> const y, MappedVector<N> ydot)>;
+    const double t, MappedConstVector<N> const& y, MappedVector<N>& ydot)>;
 
 template <unsigned N>
 using JacobianFunction = std::function<bool(const double t,
-                                            MappedConstVector<N> const y,
-                                            MappedConstVector<N> ydot,
-                                            MappedMatrix<N, N> jac)>;
+                                            MappedConstVector<N> const& y,
+                                            MappedConstVector<N> const& ydot,
+                                            MappedMatrix<N, N>& jac)>;
 
 }  // namespace MathLib
 
diff --git a/Tests/MathLib/TestCVode.cpp b/Tests/MathLib/TestCVode.cpp
index 50f54bce1217d469e494ca3cf77c74efe602912b..75b2db5a56f7aaad3eae2f5041713cd7c65f4cdb 100644
--- a/Tests/MathLib/TestCVode.cpp
+++ b/Tests/MathLib/TestCVode.cpp
@@ -16,8 +16,8 @@ const double abs_tol = 1e-8;
 const double rel_tol = 1e-8;
 
 bool f(const double,
-       MathLib::MappedConstVector<1> const y,
-       MathLib::MappedVector<1> ydot)
+       MathLib::MappedConstVector<1> const& y,
+       MathLib::MappedVector<1>& ydot)
 {
 	if (y[0] <= 0.0) return false;
 
@@ -26,9 +26,9 @@ bool f(const double,
 }
 
 bool df(const double /*t*/,
-        MathLib::MappedConstVector<1> const y,
-        MathLib::MappedConstVector<1> /*ydot*/,
-        MathLib::MappedMatrix<1, 1> jac)
+        MathLib::MappedConstVector<1> const& y,
+        MathLib::MappedConstVector<1> const& /*ydot*/,
+        MathLib::MappedMatrix<1, 1>& jac)
 {
 	if (y[0] <= 0.0) return false;
 
@@ -42,8 +42,8 @@ struct ExtraData
 };
 
 bool f_extra(const double,
-             MathLib::MappedConstVector<1> const y,
-             MathLib::MappedVector<1> ydot,
+             MathLib::MappedConstVector<1> const& y,
+             MathLib::MappedVector<1>& ydot,
              ExtraData& data)
 {
 	if (y[0] <= 0.0) return false;
@@ -153,8 +153,8 @@ TEST(MathLibCVodeTest, ExponentialExtraData)
 
 	ExtraData data;
 	auto f_lambda = [&](double t,
-	                    MathLib::MappedConstVector<1> const y,
-	                    MathLib::MappedVector<1> ydot)
+	                    MathLib::MappedConstVector<1> const& y,
+	                    MathLib::MappedVector<1>& ydot)
 	{
 		return f_extra(t, y, ydot, data);
 	};