diff --git a/examples/01_batman/main.cpp b/examples/01_batman/main.cpp
index 91bfadec8..3fd21c222 100644
--- a/examples/01_batman/main.cpp
+++ b/examples/01_batman/main.cpp
@@ -51,4 +51,4 @@ int main()
       DefaultFigure::draw_box(bi, { Color::black(), Color::red(0.4) });
   }
   #endif
-}
\ No newline at end of file
+}
diff --git a/examples/03_sivia/main.cpp b/examples/03_sivia/main.cpp
index 8ecc1a449..68c299412 100644
--- a/examples/03_sivia/main.cpp
+++ b/examples/03_sivia/main.cpp
@@ -6,6 +6,6 @@ int main()
 {
   VectorVar x(2);
   AnalyticFunction f { {x}, sqr(x[0])*sin(sqr(x[0])+sqr(x[1]))-sqr(x[1]) };
-  auto p = sivia({{-5,5},{-4,4}}, f, {0,oo}, 1e-2, true);
+  auto p = sivia_multithread({{-5,5},{-4,4}}, f, {0,oo}, 1e-2, true);
   DefaultFigure::draw_paving(p);
 }
\ No newline at end of file
diff --git a/src/core/paver/codac2_pave.cpp b/src/core/paver/codac2_pave.cpp
index 987cff977..836fcb2e5 100644
--- a/src/core/paver/codac2_pave.cpp
+++ b/src/core/paver/codac2_pave.cpp
@@ -8,6 +8,7 @@
  */
 
 #include "codac2_pave.h"
+#include <chrono>
 
 using namespace std;
 using namespace codac2;
@@ -23,10 +24,18 @@ namespace codac2
   PavingOut pave(const IntervalVector& x0, const CtcBase<IntervalVector>& c, double eps, bool verbose)
   {
     double time = 0;
-    return pave(x0,c,eps,time,verbose);
+    return pave(x0, c, eps, time, verbose);
   }
 
   PavingOut pave(const IntervalVector& x0, const CtcBase<IntervalVector>& c, double eps, double& time, bool verbose)
+  {
+    if (nb_threads()==1)
+      return pave_monothread(x0, c, eps, time, verbose);
+    else
+      return pave_multithread(x0, c, eps, time, verbose);
+  }
+
+  PavingOut pave_monothread(const IntervalVector& x0, const CtcBase<IntervalVector>& c, double eps, double& time, bool verbose)
   {
     assert_release(eps > 0.);
     assert_release(!x0.is_empty());
@@ -72,6 +81,87 @@ namespace codac2
     return p;
   }
   
+  PavingOut pave_multithread(const IntervalVector& x0, const CtcBase<IntervalVector>& c, double eps, double& time, bool verbose)
+  {
+    assert_release(eps > 0.);
+    assert_release(!x0.is_empty());
+    
+    auto start_time = std::chrono::high_resolution_clock::now();
+
+    int nthreads = nb_threads();
+
+    PavingOut p(x0);
+    // In order to be able to reconstruct the initial box, the first level represents the
+    // initial domain x0 (the left node is x0, the right one is an empty box).
+    p.tree()->bisect();
+    p.tree()->left()->boxes() = { x0 };
+    get<0>(p.tree()->right()->boxes()).set_empty();
+
+    std::shared_ptr<PavingOut_Node> n;
+    list<std::shared_ptr<PavingOut_Node>> l { p.tree()->left() };
+
+    while (l.size() < static_cast<std::size_t>(nthreads))
+    {
+      n = l.front();
+      l.pop_front();
+
+      c.contract(get<0>(n->boxes()));
+
+      if(!get<0>(n->boxes()).is_empty())
+      {
+        if(get<0>(n->boxes()).max_diam() > eps)
+        {
+          n->bisect();
+          l.push_back(n->left());
+          l.push_back(n->right());
+        }
+      }
+    }
+
+    std::vector<std::shared_ptr<PavingOut_Node>> l_vec(l.begin(), l.end());
+
+    auto worker = [&](int tid) 
+    {
+      auto xi = l_vec[tid];
+      std::list<std::shared_ptr<PavingOut_Node>> li = { xi };
+      while(!li.empty())
+      {
+        auto ni = li.front();
+        li.pop_front();
+
+        c.contract(get<0>(ni->boxes()));
+        
+        if(!get<0>(ni->boxes()).is_empty())
+        {
+          if(get<0>(ni->boxes()).max_diam() > eps)
+          {
+            ni->bisect();
+            li.push_back(ni->left());
+            li.push_back(ni->right());
+          }
+        }
+      }
+    };
+
+    std::vector<std::thread> threads;
+    for (int tid = 0; tid < nthreads; tid++)
+      threads.emplace_back(worker, tid);      
+
+    for (auto& th : threads) th.join();
+
+    std::chrono::duration<double> elapsed = std::chrono::high_resolution_clock::now() - start_time;
+
+    time = elapsed.count();
+
+    if(verbose)
+    {
+      printf("Number of thread used: %d\n", nthreads);
+      printf("Computation time: %.4fs\n\n", time);    
+    }
+
+    return p;
+  }
+
   PavingInOut pave(const IntervalVector& x0, std::shared_ptr<const SepBase> s,
     double eps, bool verbose)
   {
@@ -79,6 +169,14 @@ namespace codac2
   }
 
   PavingInOut pave(const IntervalVector& x0, const SepBase& s, double eps, bool verbose)
+  {
+    if (nb_threads()==1)
+      return pave_monothread(x0, s, eps, verbose);
+    else
+      return pave_multithread(x0, s, eps, verbose);
+  }
+
+  PavingInOut pave_monothread(const IntervalVector& x0, const SepBase& s, double eps, bool verbose)
   {
     assert_release(eps > 0.);
     assert_release(!x0.is_empty());
@@ -111,6 +209,75 @@ namespace codac2
     return p;
   }
 
+  PavingInOut pave_multithread(const IntervalVector& x0, const SepBase& s, double eps, bool verbose)
+  {
+    assert_release(eps > 0.);
+    assert_release(!x0.is_empty());
+    
+    auto start_time = std::chrono::high_resolution_clock::now();
+
+    int nthreads = nb_threads();
+
+    PavingInOut p(x0);
+    std::shared_ptr<PavingInOut_Node> n;
+    list<std::shared_ptr<PavingInOut_Node>> l { p.tree() };
+
+    while (l.size() < static_cast<std::size_t>(nthreads))
+    {
+      n = l.front();
+      l.pop_front();
+
+      auto xs = s.separate(get<0>(n->boxes()));
+      auto boundary = (xs.inner & xs.outer);
+      n->boxes() = { xs.outer, xs.inner };
+
+      if(!boundary.is_empty() && boundary.max_diam() > eps)
+      {
+        n->bisect();
+        l.push_back(n->left());
+        l.push_back(n->right());
+      }
+    }
+
+    std::vector<std::shared_ptr<PavingInOut_Node>> l_vec(l.begin(), l.end());
+
+    auto worker = [&](int tid) 
+    {
+      auto xi = l_vec[tid];
+      std::list<std::shared_ptr<PavingInOut_Node>> li = { xi };
+      while(!li.empty())
+      {
+        auto ni = li.front();
+        li.pop_front();
+
+        auto xs = s.separate(get<0>(ni->boxes()));
+        auto boundary = (xs.inner & xs.outer);
+        ni->boxes() = { xs.outer, xs.inner };
+
+        if(!boundary.is_empty() && boundary.max_diam() > eps)
+        {
+          ni->bisect();
+          li.push_back(ni->left());
+          li.push_back(ni->right());
+        }
+      }
+    };
+
+    std::vector<std::thread> threads;
+    for (int tid = 0; tid < nthreads; tid++)
+      threads.emplace_back(worker, tid);      
+
+    for (auto& th : threads) th.join();
+
+    if(verbose)
+    {
+      printf("Number of thread used: %d\n", nthreads);
+      std::chrono::duration<double> elapsed = std::chrono::high_resolution_clock::now() - start_time;
+      printf("Computation time: %.4fs\n\n", elapsed.count());    
+    }
+    return p;
+  }
+
   PavingInOut regular_pave(const IntervalVector& x0,
     const std::function<BoolInterval(const IntervalVector&)>& test,
     double eps, bool verbose)
@@ -156,6 +323,92 @@ namespace codac2
     return p;
   }
 
+  PavingInOut regular_pave_multithread(const IntervalVector& x0,
+    const std::function<BoolInterval(const IntervalVector&)>& test,
+    double eps, bool verbose)
+  {
+    assert_release(eps > 0.);
+    assert_release(!x0.is_empty());
+
+    auto start_time = std::chrono::high_resolution_clock::now();
+
+    int nthreads = nb_threads();
+
+    PavingInOut p(x0);
+    std::list<std::shared_ptr<PavingInOut_Node>> l { p.tree() };
+
+    while (l.size() < static_cast<std::size_t>(nthreads))
+    {
+      auto n = l.front();
+      l.pop_front();
+
+      auto b = test(std::get<1>(n->boxes()));
+      switch(b)
+      {
+        case BoolInterval::TRUE:
+          std::get<1>(n->boxes()).set_empty();
+          break;
+
+        case BoolInterval::FALSE:
+          std::get<0>(n->boxes()).set_empty();
+          break;
+
+        default:
+          n->bisect();
+          l.push_back(n->left());
+          l.push_back(n->right());
+      }
+    }
+
+    std::vector<std::shared_ptr<PavingInOut_Node>> l_vec(l.begin(), l.end());
+
+    auto worker = [&](int tid) 
+    {
+      auto xi = l_vec[tid];
+      std::list<std::shared_ptr<PavingInOut_Node>> li = { xi };
+      while(!li.empty())
+      {
+        auto ni = li.front();
+        li.pop_front();
+
+        auto b = test(std::get<1>(ni->boxes()));
+        switch(b)
+        {
+          case BoolInterval::TRUE:
+            std::get<1>(ni->boxes()).set_empty();
+            break;
+
+          case BoolInterval::FALSE:
+            std::get<0>(ni->boxes()).set_empty();
+            break;
+
+          default:
+            if(ni->unknown().max_diam() > eps)
+            {
+              ni->bisect();
+              li.push_back(ni->left());
+              li.push_back(ni->right());
+            }
+        }
+      }
+    };
+
+    std::vector<std::thread> threads;
+    for (int tid = 0; tid < nthreads; tid++)
+      threads.emplace_back(worker, tid);      
+
+    for (auto& th : threads) th.join();
+
+    if (verbose)
+    {
+      printf("Number of thread used: %d\n", nthreads);
+      std::chrono::duration<double> elapsed = std::chrono::high_resolution_clock::now() - start_time;
+      printf("Computation time: %.4fs\n\n", elapsed.count());
+    }
+
+    return p;
+  }
+
   PavingInOut pave_tube(const IntervalVector& x0, const SlicedTube<IntervalVector>& f, double eps, bool verbose)
   {
     return regular_pave(x0,
diff --git a/src/core/paver/codac2_pave.h b/src/core/paver/codac2_pave.h
index d1b2781a1..981d33b86 100644
--- a/src/core/paver/codac2_pave.h
+++ b/src/core/paver/codac2_pave.h
@@ -15,24 +15,50 @@
 #include "codac2_AnalyticFunction.h"
 #include "codac2_BoolInterval.h"
 #include "codac2_SlicedTube.h"
+#include "codac2_threading.h"
 
 namespace codac2
 {
   // eps: accuracy of the paving algorithm, the undefined boxes will have their max_diam <= eps
   
-  PavingOut pave(const IntervalVector& x0, std::shared_ptr<const CtcBase<IntervalVector>> c, double eps, bool verbose = false);
   PavingOut pave(const IntervalVector& x0, const CtcBase<IntervalVector>& c, double eps, double& time, bool verbose = false);
+  PavingOut pave(const IntervalVector& x0, std::shared_ptr<const CtcBase<IntervalVector>> c, double eps, bool verbose = false);
   PavingOut pave(const IntervalVector& x0, const CtcBase<IntervalVector>& c, double eps, bool verbose = false);
+  PavingOut pave_monothread(const IntervalVector& x0, const CtcBase<IntervalVector>& c, double eps, double& time, bool verbose = false);
+  PavingOut pave_multithread(const IntervalVector& x0, const CtcBase<IntervalVector>& c, double eps, double& time, bool verbose = false);
 
   PavingInOut pave(const IntervalVector& x0, std::shared_ptr<const SepBase> s, double eps, bool verbose = false);
   PavingInOut pave(const IntervalVector& x0, const SepBase& s, double eps, bool verbose = false);
+  PavingInOut pave_monothread(const IntervalVector& x0, const SepBase& s, double eps, bool verbose = false);
+  PavingInOut pave_multithread(const IntervalVector& x0, const SepBase& s, double eps, bool verbose = false);
 
   PavingInOut regular_pave(const IntervalVector& x0, const std::function<BoolInterval(const IntervalVector&)>& test, double eps, bool verbose = false);
+  PavingInOut regular_pave_multithread(const IntervalVector& x0, const std::function<BoolInterval(const IntervalVector&)>& test, double eps, bool verbose = false);
 
   template<typename Y>
   PavingInOut sivia(const IntervalVector& x0, const AnalyticFunction<Y>& f, const typename Y::Domain& y, double eps, bool verbose = false)
   {
-    return regular_pave(x0,
+    if (nb_threads()==1)
+    {
+      return regular_pave(x0,
+      [&y,&f](const IntervalVector& x)
+      {
+        auto eval = f.eval(x);
+
+        if(eval.is_subset(y))
+          return BoolInterval::TRUE;
+
+        else if(!eval.intersects(y))
+          return BoolInterval::FALSE;
+
+        else
+          return BoolInterval::UNKNOWN;
+      },
+      eps, verbose);
+    }
+    else
+    {
+      return regular_pave_multithread(x0,
       [&y,&f](const IntervalVector& x)
       {
         auto eval = f.eval(x);
@@ -47,6 +73,7 @@ namespace codac2
           return BoolInterval::UNKNOWN;
       },
       eps, verbose);
+    }
   }
 
   PavingInOut pave_tube(const IntervalVector& x0, const SlicedTube<IntervalVector>& f, double eps, bool verbose = false);