diff --git a/CMakeLists.txt b/CMakeLists.txt
index f32a570..ea599f3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -69,60 +69,5 @@ set(gcc_like $<OR:${gcc},${clang})
 if (MVPOLYNOMIAL_BUILD_TESTS)
   enable_testing()
 
-  set(boost_test_dependencies
-    algorithm
-    array
-    assert
-    bind
-    config
-    concept_check
-    container_hash
-    conversion
-    core
-    detail
-    describe
-    exception
-    function
-    functional
-    function_types
-    fusion
-    io
-    iterator
-    mp11
-    mpl
-    numeric_conversion
-    optional
-    predef
-    preprocessor
-    regex
-    range
-    smart_ptr
-    static_assert
-    throw_exception
-    tuple
-    typeof
-    type_traits
-    unordered
-    utility
-  )
-
-  foreach(deplib ${boost_test_dependencies})
-    FetchContent_Declare(
-      ${deplib}
-      GIT_REPOSITORY https://github.com/boostorg/${deplib}.git
-      GIT_TAG        boost-1.89.0
-      OVERRIDE_FIND_PACKAGE
-    )
-    FetchContent_MakeAvailable(${deplib})
-  endforeach()
-
-  FetchContent_Declare(
-    boost_unit_test
-    GIT_REPOSITORY https://github.com/boostorg/test.git
-    GIT_TAG        boost-1.89.0
-    OVERRIDE_FIND_PACKAGE
-  )
-  FetchContent_MakeAvailable(boost_unit_test)
-
   add_subdirectory(test)
 endif()
diff --git a/include/mvPolynomial/mvPolynomial.hpp b/include/mvPolynomial/mvPolynomial.hpp
index af04513..a174081 100644
--- a/include/mvPolynomial/mvPolynomial.hpp
+++ b/include/mvPolynomial/mvPolynomial.hpp
@@ -365,6 +365,8 @@ class MVPolynomial final {
     auto comparer = l.key_comp();
 
     auto mul = MVPolynomial(l.get_allocator());
+    // Clear
+    mul.erase(mul.begin());
     // Calculate all product of each l's term and r's term.
     for (const auto& l_p : l) {
       const auto& [l_idx, l_v] = l_p;
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 5af28e1..d803cd1 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,3 +1,11 @@
+FetchContent_Declare(
+  Catch2
+  GIT_REPOSITORY https://github.com/catchorg/Catch2.git
+  GIT_TAG        b3fb4b9feafcd8d91c5cb510a4775143fdbef02f # v3.11.0
+  OVERRIDE_FIND_PACKAGE
+)
+FetchContent_MakeAvailable(Catch2)
+
 add_executable(mvPolynomial_test mvPolynomial_test.cpp)
 target_compile_options(
   mvPolynomial_test
@@ -10,7 +18,7 @@ target_link_libraries(
   mvPolynomial_test
   PRIVATE
     ${PROJECT_NAME}
-    Boost::unit_test_framework
+    Catch2::Catch2WithMain
 )
 
 add_test(NAME mvPolynomial_test COMMAND mvPolynomial_test)
diff --git a/test/mvPolynomial_test.cpp b/test/mvPolynomial_test.cpp
index 5463ebb..0aa3175 100644
--- a/test/mvPolynomial_test.cpp
+++ b/test/mvPolynomial_test.cpp
@@ -1,193 +1,253 @@
-#define BOOST_TEST_MODULE mvPolynomial_unit_test
-
-#include "boost/test/unit_test.hpp"
 #include "mvPolynomial/mvPolynomial.hpp"
 
+#include <cstddef>
 #include <vector>
 
-namespace utf = boost::unit_test;
-namespace tt  = boost::test_tools;
+#include <catch2/catch_test_macros.hpp>
 
 using MP2 = mvPolynomial::MVPolynomial<int, double, 2>;
 using MP3 = mvPolynomial::MVPolynomial<int, double, 3>;
 
-BOOST_AUTO_TEST_CASE(mvPolynomial_init, *utf::tolerance(tt::fpc::percent_tolerance(1e-10))) {
-  auto ans = std::vector<std::pair<Eigen::Array2i, double>>();
-  ans      = {
-      {{0, 0}, 1},
-      {{1, 0}, 2},
-      {{0, 1}, 3},
-      {{1, 1}, 4},
-      {{2, 0}, 5},
-      {{0, 2}, 6},
-  };
+TEST_CASE("constructor", "[mvPolynomial]") {
+  SECTION("default") {
+    auto ans = std::vector<std::pair<Eigen::Array2i, double>>({
+        {{0, 0}, 0},
+    });
+    auto m   = MP2();
 
-  auto m = MP2(ans.begin(), ans.end());
-  for (auto i = 0; i < ans.size(); ++i) {
-    BOOST_TEST(m.at(ans[i].first) == ans[i].second);
+    REQUIRE(m.size() == ans.size());
+    for (auto i = 0; i < ans.size(); ++i) {
+      REQUIRE(m[ans[i].first] == ans[i].second);
+    }
+  }
+
+  SECTION("range") {
+    auto ans = std::vector<std::pair<Eigen::Array2i, double>>({
+        {{0, 0}, 1},
+        {{1, 0}, 2},
+        {{0, 1}, 3},
+        {{1, 1}, 4},
+        {{2, 0}, 5},
+        {{0, 2}, 6},
+    });
+
+    auto m = MP2(ans.begin(), ans.end());
+
+    REQUIRE(m.size() == ans.size());
+    for (size_t i = 0; i < ans.size(); ++i) {
+      REQUIRE(m.at(ans[i].first) == ans[i].second);
+    }
+  }
+
+  SECTION("initializer_list") {
+    auto ans = std::vector<std::pair<Eigen::Array2i, double>>({
+        {{0, 0}, 1},
+        {{1, 0}, 2},
+        {{0, 1}, 3},
+        {{1, 1}, 4},
+        {{2, 0}, 5},
+        {{0, 2}, 6},
+    });
+
+    auto m = MP2({
+        {{0, 0}, 1},
+        {{1, 0}, 2},
+        {{0, 1}, 3},
+        {{1, 1}, 4},
+        {{2, 0}, 5},
+        {{0, 2}, 6},
+    });
+
+    REQUIRE(m.size() == ans.size());
+    for (size_t i = 0; i < ans.size(); ++i) {
+      REQUIRE(m.at(ans[i].first) == ans[i].second);
+    }
   }
 }
 
-BOOST_AUTO_TEST_CASE(mvPolynomial_Of, *utf::tolerance(tt::fpc::percent_tolerance(1e-10))) {
-  auto ans = std::vector<std::pair<Eigen::Array2i, double>>();
-  ans      = {
+TEST_CASE("operator()", "[mvPolynomial]") {
+  auto m = MP2({
       {{0, 0}, 1},
       {{1, 0}, 2},
       {{0, 1}, 3},
       {{1, 1}, 4},
       {{2, 0}, 5},
       {{0, 2}, 6},
-  };
+  });
 
-  auto m = MP2(ans.begin(), ans.end());
-  BOOST_TEST(m(Eigen::Vector2d::Zero()) == 1);
-  BOOST_TEST(m(Eigen::Vector2d({2, 3})) == 112);
+  REQUIRE(m(Eigen::Vector2d::Zero()) == 1);
+  REQUIRE(m(Eigen::Vector2d({2, 3})) == 112);
 }
 
-BOOST_AUTO_TEST_CASE(mvPolynomial_derivative, *utf::tolerance(tt::fpc::percent_tolerance(1e-10))) {
-  auto ans = std::vector<std::pair<Eigen::Array2i, double>>();
-  ans      = {
+TEST_CASE("D", "[mvPolynomial]") {
+  auto m = MP2({
       {{0, 0}, 1},
       {{1, 0}, 2},
       {{0, 1}, 3},
       {{1, 1}, 4},
       {{2, 0}, 5},
       {{0, 2}, 6},
-  };
-  auto m = MP2(ans.begin(), ans.end());
-  ans    = {
-      {{0, 0},  2},
-      {{0, 1},  4},
-      {{1, 0}, 10}
-  };
-  auto dm0 = mvPolynomial::D(m, 0);
-  for (auto i = 0; i < ans.size(); ++i) {
-    BOOST_TEST(dm0[ans[i].first] == ans[i].second);
+  });
+
+  SECTION("D_x") {
+    auto ans = std::vector<std::pair<Eigen::Array2i, double>>({
+        {{0, 0},  2},
+        {{0, 1},  4},
+        {{1, 0}, 10}
+    });
+    auto dm0 = mvPolynomial::D(m, 0);
+
+    REQUIRE(dm0.size() == ans.size());
+    for (size_t i = 0; i < ans.size(); ++i) {
+      REQUIRE(dm0[ans[i].first] == ans[i].second);
+    }
   }
-  ans = {
-      {{0, 0},  3},
-      {{1, 0},  4},
-      {{0, 1}, 12}
-  };
-  auto dm1 = mvPolynomial::D(m, 1);
-  for (auto i = 0; i < ans.size(); ++i) {
-    BOOST_TEST(dm1[ans[i].first] == ans[i].second);
+
+  SECTION("D_y") {
+    auto ans = std::vector<std::pair<Eigen::Array2i, double>>({
+        {{0, 0},  3},
+        {{1, 0},  4},
+        {{0, 1}, 12}
+    });
+    auto dm1 = mvPolynomial::D(m, 1);
+
+    REQUIRE(dm1.size() == ans.size());
+    for (size_t i = 0; i < ans.size(); ++i) {
+      REQUIRE(dm1[ans[i].first] == ans[i].second);
+    }
   }
 }
 
-BOOST_AUTO_TEST_CASE(mvPolynomial_integral, *utf::tolerance(tt::fpc::percent_tolerance(1e-10))) {
-  auto ans = std::vector<std::pair<Eigen::Array2i, double>>();
-  ans      = {
+TEST_CASE("integral", "[mvPolynomial]") {
+  auto m = MP2({
       {{0, 0}, 1},
       {{1, 0}, 2},
       {{0, 1}, 3},
       {{1, 1}, 4},
       {{2, 0}, 5},
       {{0, 2}, 6},
-  };
-  auto m = MP2(ans.begin(), ans.end());
-  ans    = {
-      {{0, 1},   1},
-      {{1, 1},   2},
-      {{0, 2}, 1.5},
-      {{1, 2},   2},
-      {{2, 1},   5},
-      {{0, 3},   2}
-  };
-  auto sm = mvPolynomial::Integrate(m, 1);
-  for (auto i = 0; i < ans.size(); ++i) {
-    BOOST_TEST(sm[ans[i].first] == ans[i].second);
+  });
+
+  SECTION("S m dx") {
+    auto ans = std::vector<std::pair<Eigen::Array2i, double>>({
+        {{1, 0},       1},
+        {{2, 0},       1},
+        {{1, 1},       3},
+        {{2, 1},       2},
+        {{3, 0}, 5.0 / 3},
+        {{1, 2},       6}
+    });
+    auto sm  = mvPolynomial::Integrate(m, 0);
+
+    REQUIRE(sm.size() == ans.size());
+    for (size_t i = 0; i < ans.size(); ++i) {
+      REQUIRE(sm[ans[i].first] == ans[i].second);
+    }
+  }
+
+  SECTION("S m dy") {
+    auto ans = std::vector<std::pair<Eigen::Array2i, double>>({
+        {{0, 1},   1},
+        {{1, 1},   2},
+        {{0, 2}, 1.5},
+        {{1, 2},   2},
+        {{2, 1},   5},
+        {{0, 3},   2}
+    });
+    auto sm  = mvPolynomial::Integrate(m, 1);
+
+    REQUIRE(sm.size() == ans.size());
+    for (size_t i = 0; i < ans.size(); ++i) {
+      REQUIRE(sm[ans[i].first] == ans[i].second);
+    }
   }
 }
 
-BOOST_AUTO_TEST_CASE(mvPolynomial_multiply, *utf::tolerance(tt::fpc::percent_tolerance(1e-10))) {
-  auto ans = std::vector<std::pair<Eigen::Array2i, double>>();
-  ans      = {
+TEST_CASE("multiply", "[mvPolynomial]") {
+  // left hand
+  auto l = MP2({
       {{0, 0}, 1},
       {{1, 0}, 2},
       {{0, 1}, 3},
-  };
-  auto l = MP2(ans.begin(), ans.end());
-  ans    = {
-      {{2, 0}, 5},
-      {{0, 2}, 7},
-  };
-  auto m = MP2(ans.begin(), ans.end());
-  ans    = {
+  });
+
+  // right hand
+  auto r = MP2({
+      {{2, 0},  5},
+      {{0, 2},  7},
+      {{1, 1}, 11},
+  });
+
+  auto ans = std::vector<std::pair<Eigen::Array2i, double>>({
       {{2, 0},  5},
       {{3, 0}, 10},
-      {{2, 1}, 15},
+      {{2, 1}, 37}, // {{2, 1}, 15} + {{2, 1}, 22}
       {{0, 2},  7},
-      {{1, 2}, 14},
+      {{1, 2}, 47}, // {{1, 2}, 14} + {{1, 2}, 33}
       {{0, 3}, 21},
-  };
-  auto prod = l * m;
-  for (auto i = 0; i < ans.size(); ++i) {
-    BOOST_TEST(prod[ans[i].first] == ans[i].second);
-  }
-}
+      {{1, 1}, 11},
+  });
 
-BOOST_AUTO_TEST_CASE(mvPolynomial_sum, *utf::tolerance(tt::fpc::percent_tolerance(1e-10))) {
-  auto ans = std::vector<std::pair<Eigen::Array2i, double>>();
-  ans      = {
-      {{0, 0}, 1},
-      {{1, 0}, 2},
-      {{0, 1}, 3},
-  };
-  auto l = MP2(ans.begin(), ans.end());
-  ans    = {
-      {{2, 0}, 5},
-      {{0, 2}, 7},
-  };
-  auto m = MP2(ans.begin(), ans.end());
-  ans    = {
-      {{0, 0}, 1},
-      {{1, 0}, 2},
-      {{0, 1}, 3},
-      {{2, 0}, 5},
-      {{0, 2}, 7},
-  };
-  auto sum = l + m;
-  for (auto i = 0; i < ans.size(); ++i) {
-    BOOST_TEST(sum[ans[i].first] == ans[i].second);
+  auto prod = l * r;
+
+  REQUIRE(prod.size() == ans.size());
+  for (size_t i = 0; i < ans.size(); ++i) {
+    REQUIRE(prod[ans[i].first] == ans[i].second);
   }
 }
 
-BOOST_AUTO_TEST_CASE(mvPolynomial_sub, *utf::tolerance(tt::fpc::percent_tolerance(1e-10))) {
-  auto ans = std::vector<std::pair<Eigen::Array2i, double>>();
-  ans      = {
+TEST_CASE("sum", "[mvPolynomial]") {
+  auto l = MP2({
       {{0, 0}, 1},
       {{1, 0}, 2},
       {{0, 1}, 3},
-  };
-  auto l = MP2(ans.begin(), ans.end());
-  ans    = {
-      {{2, 0}, 5},
-      {{0, 2}, 7},
-  };
-  auto m = MP2(ans.begin(), ans.end());
-  ans    = {
+  });
+
+  auto r = MP2({
+      {{2, 0},  5},
+      {{0, 2},  7},
+      {{1, 1}, 11},
+  });
+
+  auto ans = std::vector<std::pair<Eigen::Array2i, double>>({
       {{0, 0},  1},
       {{1, 0},  2},
       {{0, 1},  3},
-      {{2, 0}, -5},
-      {{0, 2}, -7},
-  };
-  auto sub = l - m;
-  for (auto i = 0; i < ans.size(); ++i) {
-    BOOST_TEST(sub[ans[i].first] == ans[i].second);
+      {{2, 0},  5},
+      {{0, 2},  7},
+      {{1, 1}, 11},
+  });
+  auto sum = l + r;
+
+  REQUIRE(sum.size() == ans.size());
+  for (size_t i = 0; i < ans.size(); ++i) {
+    REQUIRE(sum[ans[i].first] == ans[i].second);
   }
 }
 
-BOOST_AUTO_TEST_CASE(
-    mvPolynomial_default_value_check, *utf::tolerance(tt::fpc::percent_tolerance(1e-10))
-) {
-  auto ans = std::vector<std::pair<Eigen::Array2i, double>>();
-  ans      = {
-      {{0, 0}, 0},
-  };
-  auto m = MP2();
+TEST_CASE("sub", "[mvPolynomial]") {
+  auto l   = MP2({
+      {{0, 0}, 1},
+      {{1, 0}, 2},
+      {{0, 1}, 3},
+  });
+  auto r   = MP2({
+      {{2, 0},  5},
+      {{0, 2},  7},
+      {{1, 1}, 11},
+  });
+  auto ans = std::vector<std::pair<Eigen::Array2i, double>>({
+      {{0, 0},   1},
+      {{1, 0},   2},
+      {{0, 1},   3},
+      {{2, 0},  -5},
+      {{0, 2},  -7},
+      {{1, 1}, -11},
+  });
+  auto sub = l - r;
+
+  REQUIRE(sub.size() == ans.size());
   for (auto i = 0; i < ans.size(); ++i) {
-    BOOST_TEST(m[ans[i].first] == ans[i].second);
+    REQUIRE(sub[ans[i].first] == ans[i].second);
   }
 }