diff --git a/mini-odeint/.gitignore b/mini-odeint/.gitignore
new file mode 100644
index 0000000..95a1449
--- /dev/null
+++ b/mini-odeint/.gitignore
@@ -0,0 +1,9 @@
+.DS_Store
+.idea
+*.log
+tmp/
+
+.cache/
+compile_commands.json
+/build/
+*.o
diff --git a/mini-odeint/CMakeLists.txt b/mini-odeint/CMakeLists.txt
new file mode 100644
index 0000000..59817e6
--- /dev/null
+++ b/mini-odeint/CMakeLists.txt
@@ -0,0 +1,28 @@
+cmake_minimum_required(VERSION 3.20)
+
+project(mini-odeint)
+
+option(FORCE_FETCH_CATCH2 "Force fetching Catch2" OFF)
+
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+find_package(Catch2 3 QUIET)
+
+if(NOT TARGET Catch2::Catch2WithMain OR FORCE_FETCH_CATCH2)
+  include(FetchContent)
+  FetchContent_Declare(
+    Catch2
+    GIT_REPOSITORY https://github.com/catchorg/Catch2.git
+    GIT_TAG v3.4.0)
+  FetchContent_MakeAvailable(Catch2)
+  list(APPEND CMAKE_MODULE_PATH ${catch2_SOURCE_DIR}/extras)
+endif()
+
+include(CTest)
+include(Catch)
+
+add_executable(tests mini-odeint-tests.cpp mini-odeint.hpp)
+target_link_libraries(tests PRIVATE Catch2::Catch2WithMain)
+
+catch_discover_tests(tests)
diff --git a/mini-odeint/mini-odeint-tests.cpp b/mini-odeint/mini-odeint-tests.cpp
new file mode 100644
index 0000000..3d37ccd
--- /dev/null
+++ b/mini-odeint/mini-odeint-tests.cpp
@@ -0,0 +1,126 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/matchers/catch_matchers_range_equals.hpp>
+
+#include "mini-odeint.hpp"
+
+#include <ranges>
+
+unsigned int Factorial(unsigned int number) {
+  return number <= 1 ? number : Factorial(number - 1) * number;
+}
+
+TEST_CASE("Factorials are computed", "[factorial]") {
+  REQUIRE(Factorial(1) == 1);
+  REQUIRE(Factorial(2) == 2);
+  REQUIRE(Factorial(3) == 6);
+  REQUIRE(Factorial(10) == 3628800);
+}
+
+TEST_CASE("Mini-odeint works", "[mini-odeint]") {
+  using namespace mini_odeint;
+
+  SECTION("float") {
+    std::vector<float> times;
+    times.reserve(1001);
+    for (int i = 0; i <= 1000; ++i) {
+      times.push_back(i / 1000.0f);
+    }
+
+    SECTION("scalar") {
+      std::vector<float> ys(times.size());
+
+      explicitRungeKutta(std::span(ys), std::span<const float>(times), 1.0f,
+                         1e-6f, [](auto y, auto t) { return -y; });
+
+      REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
+                     return std::abs(a - std::exp(-b)) < 1e-6;
+                   }));
+    }
+
+    SECTION("2 vector") {
+      std::vector<std::array<float, 2>> ys(times.size());
+
+      explicitRungeKutta(
+          std::span(ys), std::span<const float>(times), std::array{1.0f, 2.0f},
+          1e-6f, [](auto y, auto t) { return std::array{-y[0], -y[1]}; });
+
+      REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
+                     const auto s = std::exp(-b);
+
+                     return func::inf_norm(std::array{a[0] - s, a[1] - 2 * s}) <
+                            1e-6;
+                   }));
+    }
+
+    SECTION("4 vector") {
+      std::vector<std::array<float, 4>> ys(times.size());
+
+      explicitRungeKutta(std::span(ys), std::span<const float>(times),
+                         std::array{1.0f, 2.0f, 3.0f, 4.0f}, 1e-6f,
+                         [](auto y, auto t) {
+                           return std::array{-y[0], -y[1], -y[2], -y[3]};
+                         });
+
+      REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
+                     const auto s = std::exp(-b);
+
+                     return func::inf_norm(std::array{a[0] - s, a[1] - 2 * s,
+                                                      a[2] - 3 * s,
+                                                      a[3] - 4 * s}) < 1e-6;
+                   }));
+    }
+  }
+
+  SECTION("float") {
+    std::vector<double> times;
+    times.reserve(1001);
+    for (int i = 0; i <= 1000; ++i) {
+      times.push_back(i / 1000.0);
+    }
+
+    SECTION("scalar") {
+      std::vector<double> ys(times.size());
+
+      explicitRungeKutta(std::span(ys), std::span<const double>(times), 1.0,
+                         1e-6, [](auto y, auto t) { return -y; });
+
+      REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
+                     return std::abs(a - std::exp(-b)) < 1e-6;
+                   }));
+    }
+
+    SECTION("2 vector") {
+      std::vector<std::array<double, 2>> ys(times.size());
+
+      explicitRungeKutta(
+          std::span(ys), std::span<const double>(times), std::array{1.0, 2.0},
+          1e-6, [](auto y, auto t) { return std::array{-y[0], -y[1]}; });
+
+      REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
+                     const auto s = std::exp(-b);
+
+                     return func::inf_norm(std::array{a[0] - s, a[1] - 2 * s}) <
+                            1e-6;
+                   }));
+    }
+
+    SECTION("4 vector") {
+      std::vector<std::array<double, 4>> ys(times.size());
+
+      explicitRungeKutta(std::span(ys), std::span<const double>(times),
+                         std::array{1.0, 2.0, 3.0, 4.0}, 1e-6,
+                         [](auto y, auto t) {
+                           return std::array{-y[0], -y[1], -y[2], -y[3]};
+                         });
+
+      REQUIRE_THAT(ys, Catch::Matchers::RangeEquals(times, [](auto a, auto b) {
+                     const auto s = std::exp(-b);
+
+                     return func::inf_norm(std::array{a[0] - s, a[1] - 2 * s,
+                                                      a[2] - 3 * s,
+                                                      a[3] - 4 * s}) < 1e-6;
+                   }));
+    }
+  }
+}
diff --git a/mini-odeint/mini-odeint.cpp b/mini-odeint/mini-odeint.hpp
similarity index 55%
rename from mini-odeint/mini-odeint.cpp
rename to mini-odeint/mini-odeint.hpp
index 8115d4a..91d84d5 100644
--- a/mini-odeint/mini-odeint.cpp
+++ b/mini-odeint/mini-odeint.hpp
@@ -4,8 +4,6 @@
 #include <array>
 #include <cassert>
 #include <cmath>
-#include <iterator>
-#include <ranges>
 #include <span>
 
 namespace mini_odeint {
@@ -64,6 +62,127 @@ template <typename T> struct DormandPrince {
         90730570.0 / 29380423.0, -8293050.0 / 29380423.0}}};
 };
 
+template <typename T> struct scalar_type {
+  using type = T;
+};
+
+template <typename T, std::size_t N> struct scalar_type<T[N]> {
+  using type = T;
+};
+
+template <typename T, std::size_t N> struct scalar_type<std::array<T, N>> {
+  using type = T;
+};
+
+template <typename T> using scalar_type_t = typename scalar_type<T>::type;
+
+template <typename T> struct is_std_array : std::false_type {};
+template <typename T, std::size_t N>
+struct is_std_array<std::array<T, N>> : std::true_type {};
+template <typename T>
+inline constexpr bool is_std_array_v = is_std_array<T>::value;
+
+template <typename T>
+  requires std::is_trivially_copyable_v<T>
+class OdeVector {
+
+  using Scalar = scalar_type_t<T>;
+
+  T value;
+
+public:
+  OdeVector() = default;
+  explicit OdeVector(T value) : value(std::move(value)) {}
+
+  friend constexpr OdeVector operator*(OdeVector lhs, const Scalar &rhs) {
+    return lhs *= rhs;
+  }
+
+  friend constexpr OdeVector operator*(const Scalar &lhs, OdeVector rhs) {
+    return rhs *= lhs;
+  }
+
+  friend constexpr OdeVector operator+(OdeVector lhs, const OdeVector &rhs) {
+    return lhs += rhs;
+  }
+
+  friend constexpr OdeVector operator+(const T &lhs, OdeVector rhs) {
+    return rhs += OdeVector{lhs};
+  }
+
+  friend constexpr OdeVector operator+(OdeVector lhs, const T &rhs) {
+    return lhs += OdeVector{rhs};
+  }
+
+  const OdeVector &operator*=(const Scalar &rhs) {
+    if constexpr (std::is_bounded_array_v<T>) {
+      for (std::size_t i = 0; i < std::extent_v<T>; ++i) {
+        value[i] *= rhs;
+      }
+    } else if constexpr (is_std_array_v<T>) {
+      for (std::size_t i = 0; i < std::tuple_size_v<T>; ++i) {
+        value[i] *= rhs;
+      }
+    } else {
+      value *= rhs;
+    }
+    return *this;
+  }
+
+  constexpr OdeVector &operator+=(const OdeVector &rhs) {
+    if constexpr (std::is_bounded_array_v<T>) {
+      for (std::size_t i = 0; i < std::extent_v<T>; ++i) {
+        value[i] += rhs.value[i];
+      }
+    } else if constexpr (is_std_array_v<T>) {
+      for (std::size_t i = 0; i < std::tuple_size_v<T>; ++i) {
+        value[i] += rhs.value[i];
+      }
+    } else {
+      value += rhs.value;
+    }
+    return *this;
+  }
+
+  constexpr OdeVector &operator=(const T &rhs) {
+    value = rhs;
+    return *this;
+  }
+
+  explicit constexpr operator const T &() const { return value; }
+};
+
+template <typename E> struct Vec2 {
+  using value_type = E;
+
+  E x, y;
+
+  constexpr Vec2() = default;
+  constexpr Vec2(E x, E y) : x(x), y(y) {}
+
+  constexpr explicit Vec2(std::array<E, 3> v) : x(v[0]), y(v[1]) {}
+
+  friend constexpr Vec2 operator+(Vec2 lhs, const Vec2 &rhs) {
+    return lhs += rhs;
+  }
+  friend constexpr Vec2 operator*(Vec2 lhs, const value_type &rhs) {
+    return lhs *= rhs;
+  }
+  friend constexpr Vec2 operator*(const value_type &lhs, Vec2 rhs) {
+    return rhs *= lhs;
+  }
+  constexpr Vec2 &operator+=(const Vec2 &rhs) {
+    x += rhs.x;
+    y += rhs.y;
+    return *this;
+  }
+  constexpr Vec2 &operator*=(const value_type &rhs) {
+    x *= rhs;
+    y *= rhs;
+    return *this;
+  }
+};
+
 template <typename E> struct Vec3 {
   using value_type = E;
 
@@ -74,14 +193,14 @@ template <typename E> struct Vec3 {
 
   constexpr explicit Vec3(std::array<E, 3> v) : x(v[0]), y(v[1]), z(v[2]) {}
 
-  friend constexpr Vec3 operator+(const Vec3 &lhs, const Vec3 &rhs) {
-    return Vec3{lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z};
+  friend constexpr Vec3 operator+(Vec3 lhs, const Vec3 &rhs) {
+    return lhs += rhs;
   }
-  friend constexpr Vec3 operator*(const Vec3 &lhs, const value_type &rhs) {
-    return Vec3{lhs.x * rhs, lhs.y * rhs, lhs.z * rhs};
+  friend constexpr Vec3 operator*(Vec3 lhs, const value_type &rhs) {
+    return lhs *= rhs;
   }
-  friend constexpr Vec3 operator*(const value_type &lhs, const Vec3 &rhs) {
-    return Vec3{lhs * rhs.x, lhs * rhs.y, lhs * rhs.z};
+  friend constexpr Vec3 operator*(const value_type &lhs, Vec3 rhs) {
+    return rhs *= lhs;
   }
   constexpr Vec3 &operator+=(const Vec3 &rhs) {
     x += rhs.x;
@@ -89,26 +208,41 @@ template <typename E> struct Vec3 {
     z += rhs.z;
     return *this;
   }
+  constexpr Vec3 &operator*=(const value_type &rhs) {
+    x *= rhs;
+    y *= rhs;
+    z *= rhs;
+    return *this;
+  }
 };
 
-namespace alg {
-template <typename T> inline auto inf_norm(const T &v) {
-  std::ranges::max_element(v, {}, [](auto n) { return std::abs(n); });
+namespace func {
+
+inline auto inf_norm(std::ranges::range auto v) {
+  return *std::ranges::max_element(v, {}, [](auto n) { return std::abs(n); });
 }
 
-inline float inf_norm(float v) { return std::abs(v); }
+inline std::floating_point auto inf_norm(std::floating_point auto v) {
+  return std::abs(v);
+}
 
-inline double inf_norm(double v) { return std::abs(v); }
+template <typename T> inline auto inf_norm(const OdeVector<T> &v) {
+  return inf_norm(static_cast<const T &>(v));
+}
+
+template <typename E> inline E inf_norm(const Vec2<E> &v) {
+  return std::max(std::abs(v.x), std::abs(v.y));
+}
 
 template <typename E> inline E inf_norm(const Vec3<E> &v) {
   return std::max({std::abs(v.x), std::abs(v.y), std::abs(v.z)});
 }
 
-} // namespace alg
+} // namespace func
 
-template <typename Vector, typename Scalar = std::iter_value_t<Vector>,
+template <typename Vector, typename Scalar = scalar_type_t<Vector>,
           typename Tableau = DormandPrince<Scalar>>
-  requires std::same_as<Scalar, std::iter_value_t<Tableau>>
+  requires std::same_as<Scalar, typename Tableau::value_type>
 
 inline std::size_t explicitRungeKutta(std::span<Vector> ys,
                                       std::span<Scalar const> ts, Vector y0,
@@ -126,14 +260,14 @@ inline std::size_t explicitRungeKutta(std::span<Vector> ys,
   const auto &b = Tableau::b;
   const auto &b_hat = Tableau::b_hat;
 
-  static_assert(Tableau::c.back() == 1.0, "last c value must be 1.0");
+  static_assert(c.back() == 1.0, "last c value must be 1.0");
 
-  auto y_hat_n = y0;
+  OdeVector<Vector> y_hat_n{y0};
   ys[0] = y0;
 
   std::size_t it = 1;
 
-  std::array<Vector, stages> k;
+  std::array<OdeVector<Vector>, stages> k;
 
   const auto N = ts.size();
   if (!N) {
@@ -143,8 +277,6 @@ inline std::size_t explicitRungeKutta(std::span<Vector> ys,
   auto t_n = ts[0];
   auto h_n = ts[N - 1] - t_n;
 
-  int step_count = 0;
-
   k[stages - 1] = dydx(y0, t_n);
 
   while (t_n < ts[N - 1]) {
@@ -154,16 +286,17 @@ inline std::size_t explicitRungeKutta(std::span<Vector> ys,
       const auto last_k_store = k[stages - 1];
       k[0] = k[stages - 1];
       for (std::size_t i = 1; i < stages; ++i) {
-        Vector sum_ak{};
+        OdeVector<Vector> sum_ak{};
         for (std::size_t j = 0; j < i; ++j) {
           sum_ak += a[i][j] * k[j];
         }
-        k[i] = dydx(y_hat_n + h_n * sum_ak, t_n + c[i] * h_n);
+        k[i] =
+            dydx(static_cast<Vector>(y_hat_n + h_n * sum_ak), t_n + c[i] * h_n);
       }
 
       // calculate final value and error
-      Vector error{};
-      Vector sum_bk{};
+      OdeVector<Vector> error{};
+      OdeVector<Vector> sum_bk{};
       for (std::size_t i = 0; i < stages; ++i) {
         sum_bk += b_hat[i] * k[i];
         error += (b_hat[i] - b[i]) * k[i];
@@ -171,13 +304,13 @@ inline std::size_t explicitRungeKutta(std::span<Vector> ys,
       const auto y_hat_np1 = y_hat_n + h_n * sum_bk;
 
       // check if step is successful, ie error is within tolerance
-      const auto E_hp1 = alg::inf_norm(h_n * error);
+      const auto E_hp1 = func::inf_norm(h_n * error);
       if (E_hp1 < tol) {
         // if moved over any requested times then interpolate their values
         const auto t_np1 = t_n + h_n;
         while (it < N && t_np1 >= ts[it]) {
           const auto sigma = (ts[it] - t_n) / h_n;
-          Vector Phi{};
+          OdeVector<Vector> Phi{};
           for (std::size_t i = 0; i < stages; ++i) {
             auto term = sigma;
             auto b_i = term * p[i][0];
@@ -187,7 +320,7 @@ inline std::size_t explicitRungeKutta(std::span<Vector> ys,
             }
             Phi += b_i * k[i];
           }
-          ys[it] = y_hat_n + h_n * Phi;
+          ys[it] = static_cast<Vector>(y_hat_n + h_n * Phi);
           ++it;
         }
 
@@ -195,7 +328,6 @@ inline std::size_t explicitRungeKutta(std::span<Vector> ys,
         step_rejected = false;
         y_hat_n = y_hat_np1;
         t_n = t_np1;
-        ++step_count;
       } else {
         // failed step, reset last k back to stored value
         k[stages - 1] = last_k_store;
@@ -212,28 +344,4 @@ inline std::size_t explicitRungeKutta(std::span<Vector> ys,
 
 } // namespace mini_odeint
 
-#include <vector>
-
-#include <iostream>
-
-int main() {
-  using namespace mini_odeint;
-  // make vector of floats from 0.0 to 1.0 by 0.001
-  std::vector<float> times;
-  times.reserve(1001);
-  for (int i = 0; i < 1000; ++i) {
-    times.push_back(i / 1000.0);
-  }
-
-  std::vector<Vec3<float>> ys(times.size());
-
-  explicitRungeKutta(
-      std::span(ys), std::span<const float>(times), Vec3<float>{1.0, 1.0, 1.0},
-      float(1e-6), [](auto y, auto t) { return Vec3<float>{-y.x, 0.0, 0.0}; });
-
-  for (const auto &v : ys) {
-    std::cout << v.x << '\n';
-  }
-}
-
 #endif // MINI_ODEINT_H_