diff --git a/analysis/analysis.py b/analysis/analysis.py
index dd08c7d0..9fca3270 100644
--- a/analysis/analysis.py
+++ b/analysis/analysis.py
@@ -38,6 +38,8 @@ def __init__(self,label,A,b,c,b2=None):
         self._A = None
         self._P = None
         self._id = None
+        self._x_max = None
+        self._y_max = None
     
     def __iter__(self):
         """ to convert object into dictionary
@@ -155,6 +157,57 @@ def relative_error_function(self):
             self._P = sp.Abs( (self.stability_function() - sp.exp(z))/sp.exp(z) )
         return self._P
 
+    def _stability_on_real_negative_axis_explicit(self,R,x):
+        xsol = sp.reduce_abs_inequality( sp.Abs(R.evalf())-1 , '<=' , x ).as_set()
+        return sp.Abs(sp.Intersection(xsol,sp.Interval(-sp.oo,0)).start)
+
+    def _stability_on_real_negative_axis_implicit(self,R,x):
+        # force to write $R^2 - 1$ as a rational function (if implicit method)
+        expr = (1/(1/ (R**2-1) ).simplify())
+
+        if type(expr) == sp.Mul:
+            # implicit method because of rational function (sp.Mul)
+            numerator,denominator = (1,1)
+            for arg in expr.args:
+                if type(arg) == sp.Pow and arg.args[1]<0:
+                    denominator = denominator * 1/arg
+                else:
+                    numerator = numerator * arg
+        else:
+            numerator,denominator = (expr,1)
+        return sp.Abs(sp.solve_poly_inequality(sp.Poly(sp.N(numerator),x)-sp.Poly(sp.N(denominator),x),"<=")[0].evalf().start)
+
+    def stability_on_real_negative_axis(self):
+        if self._x_max is None:
+            z = sp.symbols("z")
+            x = sp.symbols("x",real=True)
+            R = self.stability_function().subs(z,x)
+            if self.is_explicit:
+                self._x_max = self._stability_on_real_negative_axis_explicit(R,x)
+            else:
+                self._x_max = self._stability_on_real_negative_axis_implicit(R,x)
+        return self._x_max
+    
+    @property
+    def x_max(self):
+        return self.stability_on_real_negative_axis()
+    
+    def stability_on_imaginary_axis(self):
+        if self._y_max is None:
+            z = sp.symbols("z")
+            y = sp.symbols("y",real=True)
+            f = sp.lambdify(y,sp.re(sp.Abs(self.stability_function().subs(z,sp.I*y))**2),'numpy')
+
+            X = np.linspace(0.,5.,1000)
+            Y = np.asarray([f(xi) for xi in X])
+            self._y_max = sp.N(min(X[Y>1.],default=sp.oo)-X[1])
+        return self._y_max
+
+    @property
+    def y_max(self):
+        return self.stability_on_imaginary_axis()
+    
+
     def scheme(self,latex=True,lawson=False):
         if latex:
             ks = sp.symbols("k_{{1:{}}}".format(self.nstages+1))
@@ -316,16 +369,39 @@ def func_name(rk_id,lawson=False):
 
         return "\n    ".join(r)
 
+def compute_box(rk):
+    def pseudo_y_max(rk):
+        z = sp.symbols("z")
+        y = sp.symbols("y",real=True)
+        expr = sp.re(sp.Abs(rk.stability_function().subs(z,sp.I*y))**2)
+        f = sp.lambdify(y,expr,'numpy')
+        X = np.linspace(0.,5.,1000)
+        Y = np.asarray([f(xi) for xi in X])
+        return sp.N(max(X[Y<=1.001],default=sp.oo))
+
+    xmax = rk.x_max
+    #use a more relaxing ymax computing
+    ymax = pseudo_y_max(rk)
+
+    maxaxis = max(xmax,ymax)
+    if maxaxis == sp.oo:
+        maxaxis = 5.0
+
+    maxaxis = 1.2*float(sp.N(maxaxis))
+    return (-maxaxis-0.5j*(maxaxis+1), 2+0.5j*(maxaxis+1))
+
 def analysis_butcher(rk):
     # to remove parenthesis around fractions in Lawson scheme expressions
     pattern = re.compile("\\\\left\(\\\\frac\{(?P<numerator>[^{}]+)\}\{(?P<denominator>[^{}]+)\}\\\\right\)")
-    
+
     analysis = {}
     # update butcher with new values
     analysis['id']          = rk.id
     analysis['nstages']     = rk.nstages
     analysis['order']       = rk.order
     analysis['stage_order'] = rk.stage_order
+    analysis['x_max']       = str(float(rk.x_max))
+    analysis['y_max']       = str(float(rk.y_max))
 
     # flags
     analysis['is_explicit'] = rk.is_explicit
@@ -392,8 +468,11 @@ def evalf_Cdomain( z , expr , zmin , zmax , N ):
 
         # add some data evaluated on complex domain
         # stability domaine
-        zmin = -5-3.5j
-        zmax =  2+3.5j
+        max_axis = max(rk.x_max,rk.y_max)
+        if max_axis == sp.oo:
+            max_axis = 5.0
+        max_axis = float(sp.N(max_axis))
+        zmin, zmax = compute_box(rk)
         N = (512,512)
         stability_domain = { 'xmin':np.real(zmin), 'xmax':np.real(zmax), 'ymin':np.imag(zmin), 'ymax':np.imag(zmax) }
         stability_domain['data'] = evalf_Cdomain( sp.symbols("z") , rk.stability_function() , zmin , zmax , N ).tolist()
diff --git a/html/filter.js b/html/filter.js
index 3fb81cce..fd147bc7 100644
--- a/html/filter.js
+++ b/html/filter.js
@@ -41,8 +41,8 @@ const sort_methods = {
   'order': (elm) => { return Number(elm.getAttribute("data-order")); } ,
   'stages': (elm) => { return Number(elm.getAttribute("data-nstages")); } ,
   'stage_order': (elm) => { return Number(elm.getAttribute("data-stage_order")); } ,
-  'xmax': (elm) => { return Number(elm.getAttribute("data-xmax")); } ,
-  'ymax': (elm) => { return Number(elm.getAttribute("data-ymax")); } ,
+  'x_max': (elm) => { return Number(elm.getAttribute("data-x_max")); } ,
+  'y_max': (elm) => { return Number(elm.getAttribute("data-y_max")); } ,
   'random': (elm) => { return Math.random() - 0.5; }
 }
 
diff --git a/html/main.js b/html/main.js
index 862d9dd3..0b5cfcc4 100644
--- a/html/main.js
+++ b/html/main.js
@@ -73,7 +73,7 @@ function stability_function(R_expr) {
   
   return div;
 }
-function resume_tableau(nstages,order,stage_order) {
+function resume_tableau(nstages,order,stage_order,x_max,y_max) {
   let tab = document.createElement("table");
   tab.classList.add("resume_tableau");
   
@@ -95,6 +95,28 @@ function resume_tableau(nstages,order,stage_order) {
   tr3.appendChild(th3); tr3.appendChild(td3);
   tab.appendChild(tr3);
 
+  let tr4 = document.createElement("tr");
+  let th4 = document.createElement("th"); th4.setAttribute('scope',"row"); th4.appendChild(document.createTextNode("stability on negative real axis"));
+  let td4 = document.createElement("td");
+  if ( x_max === Infinity ) {
+    render(String.raw`\infty`,td4);
+  } else {
+    td4.appendChild(document.createTextNode(x_max));
+  }
+  tr4.appendChild(th4); tr4.appendChild(td4);
+  tab.appendChild(tr4);
+
+  let tr5 = document.createElement("tr");
+  let th5 = document.createElement("th"); th5.setAttribute('scope',"row"); th5.appendChild(document.createTextNode("stability on imaginary axis"));
+  let td5 = document.createElement("td");
+  if ( y_max === Infinity ) {
+    render(String.raw`\infty`,td5);
+  } else {
+    td5.appendChild(document.createTextNode(y_max));
+  }
+  tr5.appendChild(th5); tr5.appendChild(td5);
+  tab.appendChild(tr5);
+
   return tab;
 }
 
@@ -366,12 +388,19 @@ function doi_bib(doi) {
 
 function rk_to_elm(rk,elm,options) {
   elm.id = rk.id;
+  if (rk.x_max === "inf") {
+    rk.x_max = Infinity;
+  }
+  if (rk.y_max === "inf") {
+    rk.y_max = Infinity;
+  }
+
   elm.classList.add('method');
   elm.setAttribute("data-order",rk.order);
   elm.setAttribute("data-nstages",rk.nstages);
   elm.setAttribute("data-stage_order",rk.stage_order);
-  elm.setAttribute("data-xmax",0.);
-  elm.setAttribute("data-ymax",0.);
+  elm.setAttribute("data-x_max",Number(rk.x_max));
+  elm.setAttribute("data-y_max",Number(rk.y_max));
   elm.setAttribute("data-explicit",rk.is_explicit);
   elm.setAttribute("data-dirk",rk.is_dirk);
   elm.setAttribute("data-embedded",rk.is_embedded);
@@ -392,7 +421,7 @@ function rk_to_elm(rk,elm,options) {
     title_R.innerHTML = "Stability function "+R.outerHTML+":";
     p.appendChild(title_R);
     p.appendChild(stability_function(rk.stability_function));
-    p.appendChild(resume_tableau(rk.nstages,rk.order,rk.stage_order));
+    p.appendChild(resume_tableau(rk.nstages,rk.order,rk.stage_order,Number(rk.x_max),Number(rk.y_max)));
     p.appendChild(prepare_canvas(rk.id));
     
     summary.addEventListener('click',function (event){
diff --git a/html/viewer.html b/html/viewer.html
index bf27e37d..84a8811d 100644
--- a/html/viewer.html
+++ b/html/viewer.html
@@ -63,8 +63,8 @@
           <option value="order" >Order</option>
           <option value="stages" >Number of stages</option>
           <option value="stage_order" >Stage order</option>
-          <option value="xmax" disabled >Real stability interval</option>
-          <option value="ymax" disabled >Imaginary stability interval</option>
+          <option value="x_max" >Real stability interval</option>
+          <option value="y_max" >Imaginary stability interval</option>
           <option value="random" >Random</option>
         </select>
         <button id="sort_btn" >Sort!</button>
diff --git a/solver/include/solver/butcher_tableau.hpp b/solver/include/solver/butcher_tableau.hpp
index 1f00f709..62bebc33 100644
--- a/solver/include/solver/butcher_tableau.hpp
+++ b/solver/include/solver/butcher_tableau.hpp
@@ -49,6 +49,9 @@ struct adaptive_butcher_tableau : public butcher_tableau<N,_value_t>
 };
 
 template <typename Tableau>
-concept is_embedded = is_butcher_tableau<Tableau> && requires (Tableau t){ t.b2; };
+concept is_embedded_tableau = is_butcher_tableau<Tableau> && requires (Tableau t){ t.b2; };
+
+template <typename T>
+concept is_embedded = is_embedded_tableau<T> || T::is_embedded == true;
 
 } // namespace ode::butcher
diff --git a/solver/include/solver/generic_butcher_rk.hpp b/solver/include/solver/generic_butcher_rk.hpp
index 97e00372..6094a4c2 100644
--- a/solver/include/solver/generic_butcher_rk.hpp
+++ b/solver/include/solver/generic_butcher_rk.hpp
@@ -4,12 +4,10 @@
 
 #include "stage.hpp"
 #include "detail.hpp"
+#include "butcher_tableau.hpp"
 
 namespace ode::butcher {
 
-template <typename Tableau>
-concept is_embedded_tableau = requires (Tableau t){ t.b2; };
-
 namespace runge_kutta {
 
   template <typename Tableau>
diff --git a/solver/include/solver/splitting.hpp b/solver/include/solver/splitting.hpp
index 387d64c3..286e974b 100644
--- a/solver/include/solver/splitting.hpp
+++ b/solver/include/solver/splitting.hpp
@@ -2,6 +2,7 @@
 
 #include <tuple>
 #include <type_traits>
+#include <concepts>
 
 namespace ode {
 namespace splitting {
@@ -30,6 +31,8 @@ namespace splitting {
   template < typename... Methods_t >
   struct lie
   {
+    static constexpr std::size_t order = 1;
+    static constexpr bool is_splitting_method = true;
     std::tuple<Methods_t...> methods;
 
     lie ( std::tuple<Methods_t...> const& t );
@@ -149,6 +152,8 @@ namespace splitting {
   template < typename... Algorithms_t >
   struct lie_tuple
   {
+    static constexpr std::size_t order = 1;
+    static constexpr bool is_splitting_method = true;
     std::tuple<Algorithms_t...> algos;
 
     lie_tuple ( Algorithms_t&&... a );
@@ -184,6 +189,8 @@ namespace splitting {
   {
     using lie<Methods_t...>::lie;
     using lie<Methods_t...>::methods;
+    static constexpr std::size_t order = 2;
+    static constexpr bool is_splitting_method = true;
 
     template < std::size_t I=0 , typename Problem_t , typename state_t , typename value_t >
     inline typename std::enable_if< (I == sizeof...(Methods_t)-1) ,void>::type
@@ -213,7 +220,6 @@ namespace splitting {
   inline typename std::enable_if< (I == sizeof...(Methods_t)-1) ,void>::type
   strang<Methods_t...>::_call_inc ( Problem_t & f , value_t tn , state_t & ui , value_t dt )
   {
-    //_call_step<I>(f,tn,ui,dt);
     ui = detail::_split_solve<I>(f,methods,ui,tn,dt);
     _call_dec<I-1>(f,tn,ui,dt);
   }
@@ -232,8 +238,7 @@ namespace splitting {
   inline typename std::enable_if< (I < sizeof...(Methods_t)-1) ,void>::type
   strang<Methods_t...>::_call_inc ( Problem_t & f , value_t tn , state_t & ui , value_t dt )
   {
-    //_call_step<I>(f,tn,ui,0.5*dt);
-    ui = detail::_split_solve<I>(f,methods,ui,tn,dt);
+    ui = detail::_split_solve<I>(f,methods,ui,tn,0.5*dt);
     _call_inc<I+1>(f,tn,ui,dt);
   }
 
@@ -243,8 +248,7 @@ namespace splitting {
   inline typename std::enable_if< (I == 0) ,void>::type
   strang<Methods_t...>::_call_dec ( Problem_t & f , value_t tn , state_t & ui , value_t dt )
   {
-    //_call_step<I>(f,tn,ui,0.5*dt);
-    ui = detail::_split_solve<I>(f,methods,ui,tn,dt);
+    ui = detail::_split_solve<I>(f,methods,ui,tn,0.5*dt);
   }
 
   /**
@@ -261,8 +265,7 @@ namespace splitting {
   inline typename std::enable_if< (I > 0) ,void>::type
   strang<Methods_t...>::_call_dec ( Problem_t & f , value_t tn , state_t & ui , value_t dt )
   {
-    //_call_step<I>(f,tn,ui,0.5*dt);
-    ui = detail::_split_solve<I>(f,methods,ui,tn,dt);
+    ui = detail::_split_solve<I>(f,methods,ui,tn,0.5*dt);
     _call_dec<I-1>(f,tn,ui,dt);
   }
 
@@ -307,6 +310,8 @@ namespace splitting {
   template < typename... Algorithms_t >
   struct strang_tuple
   {
+    static constexpr std::size_t order = 2;
+    static constexpr bool is_splitting_method = true;
     std::tuple<Algorithms_t...> algos;
 
     strang_tuple ( Algorithms_t&&... a );
@@ -332,6 +337,9 @@ namespace splitting {
     return strang_tuple<Algorithms_t...>(std::forward<Algorithms_t>(a)...);
   }
 
+  template <typename T>
+  concept is_splitting_method = requires (T t){ T::is_splitting_method == true; };
+
 
 } // namespace splitting
 } // namespace ode
diff --git a/solver/template/tpl_test_order.hxx b/solver/template/tpl_test_order.hxx
index 5ee50a81..3eb31786 100644
--- a/solver/template/tpl_test_order.hxx
+++ b/solver/template/tpl_test_order.hxx
@@ -9,12 +9,10 @@
 #include "compute_order.hpp"
 
 {% for rk in list_meth %}
-{% if rk.id != "rk_118" %}{# TODO: improve test to check order of this method #}
 TEST_CASE("order::{{ rk.id }}")
 {
   using rk_t = ode::butcher::{{ rk.id }}<>;
-  WARN( order<rk_t>() == doctest::Approx(rk_t::order).epsilon(0.05) );
+  WARN( check_order(rk_t()) == doctest::Approx(rk_t::order).epsilon(0.05) );
 }
-{% endif %}
 {% endfor %}
 
diff --git a/solver/test/compute_order.hpp b/solver/test/compute_order.hpp
index ed206dd8..d3f45da5 100644
--- a/solver/test/compute_order.hpp
+++ b/solver/test/compute_order.hpp
@@ -22,6 +22,7 @@ using namespace std::string_literals;
 
 #include <solver/solver.hpp>
 #include <solver/detail.hpp>
+#include <solver/problem.hpp>
 
 template <typename Container>
 auto
@@ -43,9 +44,23 @@ mayor_method ( Container const& x , Container const& y )
   return std::make_tuple( a, b );
 }
 
+template <typename T=double>
+auto
+error (T u, T v)
+{
+  return std::abs(u-v);
+}
+
+template <typename T=double>
+auto
+relative_error (T u, T v)
+{
+  return std::abs((u-v)/u);
+}
+
 template <typename Algorithm_t, typename T=double>
 auto
-solve_exp ( T dt , T Tf )
+solve_exp (Algorithm_t & algo, T dt , T Tf )
 {
   using state_t = T;
 
@@ -63,20 +78,13 @@ solve_exp ( T dt , T Tf )
   #else
     auto obs = [](T,state_t,T){};
   #endif
-    return ode::solve( pb , Algorithm_t(1e-5) , y0 , t_span , dt , obs );
+    return ode::solve( pb , algo , y0 , t_span , dt , obs );
 }
 
-template <typename T=double>
-auto
-error (T u, T v)
-{
-  return std::abs(u-v);
-}
 
 template <typename Algorithm_t, typename T=double>
-requires (!Algorithm_t::is_embedded)
 T
-order ()
+short_time_check_order (Algorithm_t algo=Algorithm_t())
 {
   using state_t = T;
 
@@ -93,7 +101,7 @@ order ()
 
   for ( auto n_iter : {50,25,20,15,10} ) {
     T dt = Tf/n_iter;
-    state_t u_sol = solve_exp<Algorithm_t>(dt,Tf);
+    state_t u_sol = solve_exp(algo,dt,Tf);
     auto e = error( u_exa , u_sol );
     errors.push_back(std::log( e ));
     dts.push_back(std::log(dt));
@@ -113,10 +121,68 @@ order ()
 }
 
 template <typename Algorithm_t, typename T=double>
-requires (Algorithm_t::is_embedded)
 T
-order ()
+long_time_check_order (Algorithm_t & algo)
 {
-  return Algorithm_t::order;
+  using state_t = std::valarray<T>;
+
+  T alpha=2./3., beta=4./3., gamma=1., delta=1.;
+
+  // make a problem that can be use also for splitting (in 3 parts) methods
+  auto pb = ode::make_problem(
+    [=]( T t , state_t const& u ) -> state_t {
+      return {
+        alpha*u[0] - beta*u[0]*u[1] ,
+        0.
+      };
+    },
+    [=]( T t , state_t const& u ) -> state_t {
+      return {
+        0. ,
+        delta*u[0]*u[1]
+      };
+    },
+    [=]( T t , state_t const& u ) -> state_t {
+      return {
+        0. ,
+        - gamma*u[1]
+      };
+    }
+  );
+
+  // invariant calculator
+  auto V = [=]( state_t const& u ) -> T {
+    return delta*u[0] - gamma*std::log(u[0]) + beta*u[1] - alpha*std::log(u[1]);
+  };
+
+  T x0 = 1.9;
+  state_t u_ini = {x0, x0};
+  T V_ini = V(u_ini);
+
+  std::vector<T> t_span = {0.,1000.};
+  std::vector<T> dts = {0.25,0.125,0.1,0.075,0.05}; // find a way to adapt this range to the method
+  std::vector<T> relative_errors, ldts;
+
+  for ( auto dt : dts ) {
+    state_t u_end = ode::solve(pb, algo, u_ini, t_span, dt,  [](T,state_t const&,T){});
+    relative_errors.push_back( std::log10(relative_error(V_ini,V(u_end))) );
+    ldts.push_back(std::log10(dt));
+  }
+
+  auto [a,b] = mayor_method(ldts,relative_errors);
+
+  return a;
 }
 
+template <typename Algorithm_t, typename T=double>
+T
+check_order (Algorithm_t algo=Algorithm_t())
+{
+  if constexpr ( ode::butcher::is_embedded<Algorithm_t> ) {
+    return Algorithm_t::order;
+  } else if constexpr ( Algorithm_t::order >= 8 || ode::splitting::is_splitting_method<Algorithm_t> ) {
+    return long_time_check_order(algo);
+  } else {
+    return short_time_check_order(algo);
+  }
+}
diff --git a/solver/test/main.cpp b/solver/test/main.cpp
index bcf70c68..5730e707 100644
--- a/solver/test/main.cpp
+++ b/solver/test/main.cpp
@@ -3,4 +3,5 @@
 
 #include "detail.hxx"
 #include "test_order.hxx"
+#include "test_order_split.hxx"
 
diff --git a/solver/test/test_order_split.hxx b/solver/test/test_order_split.hxx
new file mode 100644
index 00000000..1c9cd060
--- /dev/null
+++ b/solver/test/test_order_split.hxx
@@ -0,0 +1,30 @@
+#pragma once
+
+#include <vector>
+#include <cmath>
+#include <algorithm>
+#include <iostream>
+
+#include <doctest/doctest.h>
+
+#include <solver/splitting.hpp>
+#include <solver/butcher_methods.hpp>
+
+#include "compute_order.hpp"
+
+TEST_CASE("detail::power")
+{
+  auto lie_splitting = ode::splitting::make_lie_tuple(
+    ode::butcher::rk_33<>(),
+    ode::butcher::rk_33<>(),
+    ode::butcher::rk_33<>()
+  );
+  auto strang_splitting = ode::splitting::make_strang_tuple(
+    ode::butcher::rk_33<>(),
+    ode::butcher::rk_33<>(),
+    ode::butcher::rk_33<>()
+  );
+
+  WARN( check_order(lie_splitting) == doctest::Approx(lie_splitting.order).epsilon(0.05) );
+  WARN( check_order(strang_splitting) == doctest::Approx(strang_splitting.order).epsilon(0.05) );
+}