diff --git a/doc/index.rst b/doc/index.rst
index 05380335..6e4a0a48 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -360,6 +360,10 @@ The following global options are defined. Note that all are typically formatted
      - :python:`["oo", "zoo", "nan", "NaN", "__h"]`
      - list of strings
      - For each forbidden name: emit an error if a variable or parameter by this name occurs in the input.
+   * - ``use_alternative_expM``
+     - :python:`False`
+     - boolean
+     - If :python:`False`, use the sympy function ``sympy.exp`` to compute the matrix exponential. If :python:`True`, use an alternative function (see :py:func:`odetoolbox.sympy_helpers.expMt` for details). This can be useful as calls to ``sympy.exp`` can sometimes take a very large amount of time.
 
 
 Output
diff --git a/odetoolbox/__init__.py b/odetoolbox/__init__.py
index 89c42cc6..e1790245 100644
--- a/odetoolbox/__init__.py
+++ b/odetoolbox/__init__.py
@@ -36,8 +36,8 @@
     import pygsl.odeiv as odeiv
     PYGSL_AVAILABLE = True
 except ImportError as ie:
-    logging.warning("PyGSL is not available. The stiffness test will be skipped.")
-    logging.warning("Error when importing: " + str(ie))
+    logging.getLogger(__name__).warning("PyGSL is not available. The stiffness test will be skipped.")
+    logging.getLogger(__name__).warning("Error when importing: " + str(ie))
     PYGSL_AVAILABLE = False
 
 if PYGSL_AVAILABLE:
@@ -45,6 +45,7 @@
 
 try:
     import graphviz
+    logging.getLogger("graphviz").setLevel(logging.ERROR)
     PLOT_DEPENDENCY_GRAPH = True
 except ImportError:
     PLOT_DEPENDENCY_GRAPH = False
@@ -61,7 +62,7 @@ def _find_analytically_solvable_equations(shape_sys, shapes, parameters=None):
 
     Perform dependency analysis and plot dependency graph.
     """
-    logging.info("Finding analytically solvable equations...")
+    logging.getLogger(__name__).debug("Finding analytically solvable equations...")
     dependency_edges = shape_sys.get_dependency_edges()
 
     if PLOT_DEPENDENCY_GRAPH:
@@ -94,7 +95,7 @@ def _read_global_config(indict):
     r"""
     Process global configuration options.
     """
-    logging.info("Processing global options...")
+    logging.getLogger(__name__).debug("Processing global options...")
     if "options" in indict.keys():
         for key, value in indict["options"].items():
             assert key in Config.config.keys(), "Unknown key specified in global options dictionary: \"" + str(key) + "\""
@@ -108,7 +109,7 @@ def _from_json_to_shapes(indict, parameters=None) -> Tuple[List[Shape], Dict[sym
     :param indict: ODE-toolbox input dictionary.
     """
 
-    logging.info("Processing input shapes...")
+    logging.getLogger(__name__).debug("Processing input...")
 
     # first run for grabbing all the variable names. Coefficients might be incorrect.
     all_variable_symbols = []
@@ -122,7 +123,6 @@ def _from_json_to_shapes(indict, parameters=None) -> Tuple[List[Shape], Dict[sym
     all_parameter_symbols -= all_variable_symbols_
     del all_variable_symbols_
     assert all([_is_sympy_type(sym) for sym in all_variable_symbols])
-    logging.info("All known variables: " + str(all_variable_symbols) + ", all parameters used in ODEs: " + str(all_parameter_symbols))
 
     # validate input for forbidden names
     for var in set(all_variable_symbols) | all_parameter_symbols:
@@ -136,13 +136,13 @@ def _from_json_to_shapes(indict, parameters=None) -> Tuple[List[Shape], Dict[sym
         assert isinstance(param, SympyExpr)
         if not param in parameters.keys():
             # this parameter was used in an ODE, but not explicitly numerically specified
-            logging.info("No numerical value specified for parameter \"" + str(param) + "\"")    # INFO level because this is OK!
+            logging.getLogger(__name__).info("No numerical value specified for parameter \"" + str(param) + "\"")    # INFO level because this is OK!
             parameters[param] = None
 
     # second run with the now-known list of variable symbols
     shapes = []
     for shape_json in indict["dynamics"]:
-        shape = Shape.from_json(shape_json, all_variable_symbols=all_variable_symbols, parameters=parameters, _debug=True)
+        shape = Shape.from_json(shape_json, all_variable_symbols=all_variable_symbols, parameters=parameters)
         shapes.append(shape)
 
     return shapes, parameters
@@ -193,11 +193,11 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so
 
     _init_logging(log_level)
 
-    logging.info("Analysing input:")
-    logging.info(json.dumps(indict, indent=4, sort_keys=True))
+    logging.getLogger(__name__).info("Analysing input:")
+    logging.getLogger(__name__).info(json.dumps(indict, indent=4, sort_keys=True))
 
     if "dynamics" not in indict:
-        logging.info("Warning: empty input (no dynamical equations found); returning empty output")
+        logging.getLogger(__name__).info("Warning: empty input (no dynamical equations found); returning empty output")
         return [], SystemOfShapes.from_shapes([]), []
 
     _read_global_config(indict)
@@ -224,17 +224,17 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so
 
     for shape in shapes:
         if not shape.is_homogeneous() and shape.order > 1:
-            logging.error("For symbol " + str(shape.symbol) + ": higher-order inhomogeneous ODEs are not supported")
+            logging.getLogger(__name__).error("For symbol " + str(shape.symbol) + ": higher-order inhomogeneous ODEs are not supported")
             sys.exit(1)
 
     shape_sys = SystemOfShapes.from_shapes(shapes, parameters=parameters)
     _, node_is_analytically_solvable = _find_analytically_solvable_equations(shape_sys, shapes, parameters=parameters)
 
-    logging.debug("System of equations:")
-    logging.debug("x = " + str(shape_sys.x_))
-    logging.debug("A = " + repr(shape_sys.A_))
-    logging.debug("b = " + str(shape_sys.b_))
-    logging.debug("c = " + str(shape_sys.c_))
+    logging.getLogger(__name__).info("System of equations (with dx/dt = Ax + b + c):")
+    logging.getLogger(__name__).info("x = " + str(shape_sys.x_))
+    logging.getLogger(__name__).info("A = " + repr(shape_sys.A_))
+    logging.getLogger(__name__).info("b = " + str(shape_sys.b_))
+    logging.getLogger(__name__).info("c = " + str(shape_sys.c_))
 
 
     #
@@ -249,7 +249,7 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so
 
     analytic_solver_json = None
     if analytic_syms:
-        logging.info("Generating propagators for the following symbols: " + ", ".join([str(k) for k in analytic_syms]))
+        logging.getLogger(__name__).info("Generating propagators for the following symbols: " + ", ".join([str(k) for k in analytic_syms]))
         sub_sys = shape_sys.get_sub_system(analytic_syms)
         analytic_solver_json = sub_sys.generate_propagator_solver(disable_singularity_detection=disable_singularity_detection)
         analytic_solver_json["solver"] = "analytical"
@@ -262,7 +262,7 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so
 
     if len(analytic_syms) < len(shape_sys.x_):
         numeric_syms = list(set(shape_sys.x_) - set(analytic_syms))
-        logging.info("Generating numerical solver for the following symbols: " + ", ".join([str(sym) for sym in numeric_syms]))
+        logging.getLogger(__name__).info("Generating numerical solver for the following symbols: " + ", ".join([str(sym) for sym in numeric_syms]))
         sub_sys = shape_sys.get_sub_system(numeric_syms)
         solver_json = sub_sys.generate_numeric_solver(state_variables=shape_sys.x_)
         solver_json["solver"] = "numeric"   # will be appended to if stiffness testing is used
@@ -270,7 +270,7 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so
             if not PYGSL_AVAILABLE:
                 raise Exception("Stiffness test requested, but PyGSL not available")
 
-            logging.info("Performing stiffness test...")
+            logging.getLogger(__name__).info("Performing stiffness test...")
             kwargs = {}   # type: Dict[str, Any]
             if "options" in indict.keys() and "random_seed" in indict["options"].keys():
                 random_seed = int(indict["options"]["random_seed"])
@@ -289,7 +289,7 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so
             solver_type = tester.check_stiffness()
             if not solver_type is None:
                 solver_json["solver"] += "-" + solver_type
-                logging.info(solver_type + " scheme")
+                logging.getLogger(__name__).info(solver_type + " scheme")
 
         solvers_json.append(solver_json)
 
@@ -376,10 +376,10 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so
 
                 if preserve_expressions and sym in preserve_expressions:
                     if "analytic" in solver_json["solver"]:
-                        logging.warning("Not preserving expression for variable \"" + sym + "\" as it is solved by propagator solver")
+                        logging.getLogger(__name__).warning("Not preserving expression for variable \"" + sym + "\" as it is solved by propagator solver")
                         continue
 
-                    logging.info("Preserving expression for variable \"" + sym + "\"")
+                    logging.getLogger(__name__).info("Preserving expression for variable \"" + sym + "\"")
                     var_def_str = _find_variable_definition(indict, sym, order=1)
                     assert var_def_str is not None
                     solver_json["update_expressions"][sym] = var_def_str.replace("'", Config().differential_order_symbol)
@@ -388,8 +388,8 @@ def _analysis(indict, disable_stiffness_check: bool = False, disable_analytic_so
             for sym, expr in solver_json["propagators"].items():
                 solver_json["propagators"][sym] = str(expr)
 
-    logging.info("In ode-toolbox: returning outdict = ")
-    logging.info(json.dumps(solvers_json, indent=4, sort_keys=True))
+    logging.getLogger(__name__).info("Final output result:")
+    logging.getLogger(__name__).info(json.dumps(solvers_json, indent=4, sort_keys=True))
 
     return solvers_json, shape_sys, shapes
 
@@ -400,9 +400,10 @@ def _init_logging(log_level: Union[str, int] = logging.WARNING):
 
     :param log_level: Sets the logging threshold. Logging messages which are less severe than ``log_level`` will be ignored. Log levels can be provided as an integer or string, for example "INFO" (more messages) or "WARN" (fewer messages). For a list of valid logging levels, see https://docs.python.org/3/library/logging.html#logging-levels
     """
-    fmt = '%(levelname)s:%(message)s'
+    fmt = "[ODE-toolbox] %(levelname)s:%(message)s"
     logging.basicConfig(format=fmt)
-    logging.getLogger().setLevel(log_level)
+    logger = logging.getLogger(__name__)
+    logger.setLevel(log_level)
 
 
 def analysis(indict, disable_stiffness_check: bool = False, disable_analytic_solver: bool = False, disable_singularity_detection: bool = False, preserve_expressions: Union[bool, Iterable[str]] = False, log_level: Union[str, int] = logging.WARNING) -> List[Dict]:
diff --git a/odetoolbox/config.py b/odetoolbox/config.py
index 09444a6c..43f3f6a5 100644
--- a/odetoolbox/config.py
+++ b/odetoolbox/config.py
@@ -36,7 +36,8 @@ class Config:
         "max_step_size": 999.,
         "integration_accuracy_abs": 1E-6,
         "integration_accuracy_rel": 1E-6,
-        "forbidden_names": ["oo", "zoo", "nan", "NaN", "__h"]
+        "forbidden_names": ["oo", "zoo", "nan", "NaN", "__h"],
+        "use_alternative_expM": False
     }
 
     def __getitem__(self, key):
diff --git a/odetoolbox/dependency_graph_plotter.py b/odetoolbox/dependency_graph_plotter.py
index 0256d1c0..546218af 100644
--- a/odetoolbox/dependency_graph_plotter.py
+++ b/odetoolbox/dependency_graph_plotter.py
@@ -75,5 +75,5 @@ def plot_graph(cls, shapes, dependency_edges, node_is_lin, fn=None):
             dot.edge(str(e[0]), str(e[1]))
 
         if not fn is None:
-            logging.info("Saving dependency graph plot to " + fn)
+            logging.getLogger(__name__).debug("Saving dependency graph plot to " + fn)
             dot.render(fn)
diff --git a/odetoolbox/mixed_integrator.py b/odetoolbox/mixed_integrator.py
index 9509d7f6..7a4c134c 100644
--- a/odetoolbox/mixed_integrator.py
+++ b/odetoolbox/mixed_integrator.py
@@ -42,8 +42,8 @@
     import pygsl.odeiv as odeiv
     PYGSL_AVAILABLE = True
 except ImportError as ie:
-    logging.warning("PyGSL is not available. The stiffness test will be skipped.")
-    logging.warning("Error when importing: " + str(ie))
+    logging.getLogger(__name__).warning("PyGSL is not available. The stiffness test will be skipped.")
+    logging.getLogger(__name__).warning("Error when importing: " + str(ie))
     PYGSL_AVAILABLE = False
 
 
@@ -268,7 +268,7 @@ def integrate_ode(self, initial_values=None, h_min_lower_bound=5E-9, raise_error
 
                     if h_min < h_min_lower_bound:
                         estr = "Integration step below %.e (s=%.f). Please check your ODE." % (h_min_lower_bound, h_min)
-                        logging.warning(estr)
+                        logging.getLogger(__name__).warning(estr)
                         if raise_errors:
                             raise Exception(estr)
 
@@ -345,7 +345,7 @@ def integrate_ode(self, initial_values=None, h_min_lower_bound=5E-9, raise_error
             if self._debug_plot_dir:
                 self.integrator_debug_plot(t_log, h_log, y_log, dir=self._debug_plot_dir)
 
-        logging.info("For integrator = " + str(self.numeric_integrator) + ": h_min = " + str(h_min) + ", h_avg = " + str(h_avg) + ", runtime = " + str(runtime))
+        logging.getLogger(__name__).info("For integrator = " + str(self.numeric_integrator) + ": h_min = " + str(h_min) + ", h_avg = " + str(h_avg) + ", runtime = " + str(runtime))
 
         sym_list = self._system_of_shapes.x_
 
@@ -471,10 +471,10 @@ def step(self, t, y, params):
             # return [ float(self._update_expr[str(sym)].evalf(subs=self._locals)) for sym in self._system_of_shapes.x_ ]	# non-wrapped version
             _ret = [self._update_expr_wrapped[str(sym)](*y) for sym in self._system_of_shapes.x_]
         except Exception as e:
-            logging.error("E==>", type(e).__name__ + ": " + str(e))
-            logging.error("     Local parameters at time of failure:")
+            logging.getLogger(__name__).error("E==>", type(e).__name__ + ": " + str(e))
+            logging.getLogger(__name__).error("     Local parameters at time of failure:")
             for k, v in self._locals.items():
-                logging.error("    ", k, "=", v)
+                logging.getLogger(__name__).error("    ", k, "=", v)
             raise
 
         return _ret
diff --git a/odetoolbox/shapes.py b/odetoolbox/shapes.py
index 0cbf3ef6..97ec6c9d 100644
--- a/odetoolbox/shapes.py
+++ b/odetoolbox/shapes.py
@@ -160,8 +160,6 @@ def __init__(self, symbol, order, initial_values, derivative_factors, inhom_term
         if not self.upper_bound is None:
             self.upper_bound = _custom_simplify_expr(self.upper_bound)
 
-        logging.debug("Created Shape with symbol " + str(self.symbol) + ", derivative_factors = " + str(self.derivative_factors) + ", inhom_term = " + str(self.inhom_term) + ", nonlin_term = " + str(self.nonlin_term))
-
 
     def __str__(self):
         s = "Shape \"" + str(self.symbol) + "\" of order " + str(self.order)
@@ -267,7 +265,7 @@ def _parse_defining_expression(cls, s: str) -> Tuple[str, int, str]:
 
 
     @classmethod
-    def from_json(cls, indict, all_variable_symbols=None, parameters=None, _debug=False):
+    def from_json(cls, indict, all_variable_symbols=None, parameters=None):
         r"""
         Create a :python:`Shape` instance from an input dictionary.
 
@@ -346,7 +344,7 @@ def reconstitute_expr(self) -> SympyExpr:
         derivative_symbols = self.get_state_variables(derivative_symbol=Config().differential_order_symbol)
         for derivative_factor, derivative_symbol in zip(self.derivative_factors, derivative_symbols):
             expr += derivative_factor * derivative_symbol
-        logging.debug("Shape " + str(self.symbol) + ": reconstituting expression " + str(expr))
+
         return expr
 
 
@@ -371,8 +369,6 @@ def split_lin_inhom_nonlin(expr, x, parameters=None):
         if parameters is None:
             parameters = {}
 
-        logging.debug("Splitting expression " + str(expr) + " (symbols " + str(x) + ")")
-
         lin_factors = sympy.zeros(len(x), 1)
         inhom_term = sympy.Float(0)
         nonlin_term = sympy.Float(0)
@@ -397,15 +393,11 @@ def split_lin_inhom_nonlin(expr, x, parameters=None):
                 if not is_lin:
                     nonlin_term += term
 
-        logging.debug("\tlinear factors: " + str(lin_factors))
-        logging.debug("\tinhomogeneous term: " + str(inhom_term))
-        logging.debug("\tnonlinear term: " + str(nonlin_term))
-
         return lin_factors, inhom_term, nonlin_term
 
 
     @classmethod
-    def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_variable_symbols=None, debug=False) -> Shape:
+    def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_variable_symbols=None) -> Shape:
         r"""
         Create a Shape object given a function of time.
 
@@ -434,9 +426,6 @@ def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_vari
         # `derivatives` is a list of all derivatives of `shape` up to the order we are checking, starting at 0.
         derivatives = [definition, sympy.diff(definition, Config().input_time_symbol)]
 
-        logging.debug("Processing function-of-time shape \"" + symbol + "\" with defining expression = \"" + str(definition) + "\"")
-
-
         #
         #   to avoid a division by zero below, we have to find a `t` so that the shape function is not zero at this `t`.
         #
@@ -447,8 +436,6 @@ def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_vari
                 t_val = t_
                 break
 
-        logging.debug("Found t: " + str(t_val))
-
         if t_val is None:
 
             #
@@ -465,8 +452,6 @@ def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_vari
 
         order = 1
 
-        logging.debug("\tFinding ode for order 1...")
-
         derivative_factors = [(1 / derivatives[0] * derivatives[1]).subs(Config().input_time_symbol, t_val)]
         diff_rhs_lhs = derivatives[1] - derivative_factors[0] * derivatives[0]
         found_ode = _is_zero(diff_rhs_lhs)
@@ -479,8 +464,6 @@ def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_vari
         while not found_ode and order < max_order:
             order += 1
 
-            logging.debug("\tFinding ode for order " + str(order) + "...")
-
             # Add the next higher derivative to the list
             derivatives.append(sympy.diff(derivatives[-1], Config().input_time_symbol))
 
@@ -522,7 +505,6 @@ def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_vari
             #
 
             diff_rhs_lhs = 0
-            logging.debug("\tchecking whether shape definition is satisfied...")
             for k in range(order):
                 diff_rhs_lhs -= derivative_factors[k] * derivatives[k]
             diff_rhs_lhs += derivatives[order]
@@ -534,7 +516,6 @@ def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_vari
         if not found_ode:
             raise Exception("Shape does not satisfy any ODE of order <= " + str(max_order))
 
-        logging.debug("Shape satisfies ODE of order = " + str(order))
 
         #
         #    calculate the initial values of the found ODE
@@ -546,7 +527,7 @@ def from_function(cls, symbol: str, definition, max_t=100, max_order=4, all_vari
 
 
     @classmethod
-    def from_ode(cls, symbol: str, definition: str, initial_values: dict, all_variable_symbols=None, lower_bound=None, upper_bound=None, parameters=None, debug=False, **kwargs) -> Shape:
+    def from_ode(cls, symbol: str, definition: str, initial_values: dict, all_variable_symbols=None, lower_bound=None, upper_bound=None, parameters=None, **kwargs) -> Shape:
         r"""
         Create a :python:`Shape` object given an ODE and initial values.
 
@@ -568,8 +549,6 @@ def from_ode(cls, symbol: str, definition: str, initial_values: dict, all_variab
         assert type(definition) is str
         assert type(initial_values) is dict
 
-        logging.debug("\nProcessing differential-equation form shape " + str(symbol) + " with defining expression = \"" + str(definition) + "\"")
-
         if all_variable_symbols is None:
             all_variable_symbols = []
 
@@ -605,6 +584,5 @@ def from_ode(cls, symbol: str, definition: str, initial_values: dict, all_variab
             nonlin_term = nonlin_term + functools.reduce(lambda x, y: x + y, nonlocal_derivative_terms)
 
         shape = cls(sympy.Symbol(symbol), order, initial_values, local_derivative_factors, inhom_term, nonlin_term, lower_bound, upper_bound)
-        logging.debug("\tReturning shape: " + str(shape))
 
         return shape
diff --git a/odetoolbox/singularity_detection.py b/odetoolbox/singularity_detection.py
index 4915776e..df816ddf 100644
--- a/odetoolbox/singularity_detection.py
+++ b/odetoolbox/singularity_detection.py
@@ -169,7 +169,7 @@ def find_inhomogeneous_singularities(expr) -> Set[SymmetricEq]:
         conditions
             a set with equations, where the left-hand side of each equation is the variable that is to be subsituted, and the right-hand side is the expression to put in its place
         """
-        logging.debug("Checking for singularities (divisions by zero) in the inhomogeneous part of the update equations...")
+        logging.getLogger(__name__).debug("Checking for singularities (divisions by zero) in the inhomogeneous part of the update equations...")
 
         symbols = list(expr.free_symbols)
         conditions = set()
@@ -178,10 +178,10 @@ def find_inhomogeneous_singularities(expr) -> Set[SymmetricEq]:
             if conditions:
                 # if there is one or more condition under which the solution goes to infinity...
 
-                logging.warning("Under certain conditions, one or more inhomogeneous term(s) in the system contain a division by zero.")
-                logging.warning("List of all conditions that result in a division by zero:")
+                logging.getLogger(__name__).warning("Under certain conditions, one or more inhomogeneous term(s) in the system contain a division by zero.")
+                logging.getLogger(__name__).warning("List of all conditions that result in a division by zero:")
                 for cond_set in conditions:
-                    logging.warning("\t" + r" ∧ ".join([str(eq.lhs) + " = " + str(eq.rhs) for eq in cond_set]))
+                    logging.getLogger(__name__).warning("\t" + r" ∧ ".join([str(eq.lhs) + " = " + str(eq.rhs) for eq in cond_set]))
 
         return conditions
 
@@ -203,7 +203,7 @@ def find_propagator_singularities(P: sympy.Matrix, A: sympy.Matrix) -> Set[Symme
         conditions
             a set with equations, where the left-hand side of each equation is the variable that is to be subsituted, and the right-hand side is the expression to put in its place
         """
-        logging.debug("Checking for singularities (divisions by zero) in the propagator matrix...")
+        logging.getLogger(__name__).debug("Checking for singularities (divisions by zero) in the propagator matrix...")
         try:
             conditions = SingularityDetection._generate_singularity_conditions(P)
             conditions = SingularityDetection._filter_valid_conditions(conditions, A)  # filters out the invalid conditions (invalid means those for which A is not defined)
diff --git a/odetoolbox/stiffness.py b/odetoolbox/stiffness.py
index 7de05369..7ba0f90b 100644
--- a/odetoolbox/stiffness.py
+++ b/odetoolbox/stiffness.py
@@ -34,8 +34,8 @@
     import pygsl.odeiv as odeiv
     PYGSL_AVAILABLE = True
 except ImportError as ie:
-    logging.warning("PyGSL is not available. The stiffness test will be skipped.")
-    logging.warning("Error when importing: " + str(ie))
+    logging.getLogger(__name__).warning("PyGSL is not available. The stiffness test will be skipped.")
+    logging.getLogger(__name__).warning("Error when importing: " + str(ie))
     PYGSL_AVAILABLE = False
 
 
@@ -111,10 +111,10 @@ def check_stiffness(self, raise_errors=False):
             step_min_exp, step_average_exp, runtime_exp = self._evaluate_integrator(odeiv.step_rk4, raise_errors=raise_errors)
             step_min_imp, step_average_imp, runtime_imp = self._evaluate_integrator(odeiv.step_bsimp, raise_errors=raise_errors)
         except ParametersIncompleteException:
-            logging.warning("Stiffness test not possible because numerical values were not specified for all parameters.")
+            logging.getLogger(__name__).warning("Stiffness test not possible because numerical values were not specified for all parameters.")
             return None
 
-        # logging.info("runtime (imp:exp): %f:%f" % (runtime_imp, runtime_exp))
+        # logging.getLogger(__name__).info("runtime (imp:exp): %f:%f" % (runtime_imp, runtime_exp))
 
         return self._draw_decision(step_min_imp, step_min_exp, step_average_imp, step_average_exp)
 
@@ -144,7 +144,7 @@ def _evaluate_integrator(self, integrator, h_min_lower_bound=1E-12, raise_errors
         #  initialise and run mixed integrator
         #
 
-        logging.info("Simulating for " + str(self.sim_time) + " with max_step_size = " + str(self.max_step_size))
+        logging.getLogger(__name__).info("Simulating for " + str(self.sim_time) + " with max_step_size = " + str(self.max_step_size))
 
         mixed_integrator = MixedIntegrator(integrator,
                                            self.system_of_shapes,
@@ -160,7 +160,7 @@ def _evaluate_integrator(self, integrator, h_min_lower_bound=1E-12, raise_errors
                                            alias_spikes=self.alias_spikes)
         h_min, h_avg, runtime = (lambda x: x[:3])(mixed_integrator.integrate_ode(h_min_lower_bound=h_min_lower_bound, raise_errors=raise_errors, debug=debug))
 
-        logging.info("For integrator = " + str(integrator) + ": h_min = " + str(h_min) + ", h_avg = " + str(h_avg) + ", runtime = " + str(runtime))
+        logging.getLogger(__name__).info("For integrator = " + str(integrator) + ": h_min = " + str(h_min) + ", h_avg = " + str(h_avg) + ", runtime = " + str(runtime))
 
         return h_min, h_avg, runtime
 
diff --git a/odetoolbox/sympy_helpers.py b/odetoolbox/sympy_helpers.py
index 9c1ac74c..307b0399 100644
--- a/odetoolbox/sympy_helpers.py
+++ b/odetoolbox/sympy_helpers.py
@@ -105,7 +105,7 @@ def _custom_simplify_expr(expr: str):
     try:
         # skip simplification for long expressions
         if len(str(expr)) > Config().expression_simplification_threshold:
-            logging.warning("Length of expression \"" + str(expr) + "\" exceeds sympy simplification threshold")
+            logging.getLogger(__name__).warning("Length of expression \"" + str(expr) + "\" exceeds sympy simplification threshold")
 
         _simplify_expr = compile(Config().simplify_expression, filename="<string>", mode="eval")
         expr_simplified = eval(_simplify_expr)
@@ -166,3 +166,73 @@ def _print_Function(self, expr):
             return expr.func.__name__.lower() + "(%s)" % self.stringify(expr.args, ", ")
 
         return expr.func.__name__ + "(%s)" % self.stringify(expr.args, ", ")
+
+
+
+def expMt(M, t=1):
+    """Compute matrix exponential exp(M*t).
+
+    Based on code contributed by GitHub user @oscarbenjamin, July 29th, 2021 [1]_.
+
+    .. [1] https://github.com/sympy/sympy/issues/21585
+    """
+    def ilt(e, s, t):
+        """Fast inverse Laplace transform of rational function including RootSum"""
+        a, b, n = sympy.symbols('a, b, n', cls=sympy.Wild, exclude=[s])
+
+        def _ilt(e):
+            if not e.has(s):
+                return e
+            elif e.is_Add:
+                return _ilt_add(e)
+            elif e.is_Mul:
+                return _ilt_mul(e)
+            elif e.is_Pow:
+                return _ilt_pow(e)
+            elif isinstance(e, sympy.RootSum):
+                return _ilt_rootsum(e)
+            else:
+                raise NotImplementedError
+
+        def _ilt_add(e):
+            return e.func(*map(_ilt, e.args))
+
+        def _ilt_mul(e):
+            coeff, expr = e.as_independent(s)
+            if expr.is_Mul:
+                raise NotImplementedError
+            return coeff * _ilt(expr)
+
+        def _ilt_pow(e):
+            match = e.match((a*s + b)**n)
+            if match is not None:
+                nm, am, bm = match[n], match[a], match[b]
+                if nm.is_Integer and nm < 0:
+                    if nm == 1:
+                        return sympy.exp(-(bm/am)*t) / am
+                    else:
+                        return t**(-nm-1)*sympy.exp(-(bm/am)*t)/(am**-nm*sympy.gamma(-nm))
+            raise NotImplementedError
+
+        def _ilt_rootsum(e):
+            expr = e.fun.expr
+            [variable] = e.fun.variables
+            return sympy.RootSum(e.poly, sympy.Lambda(variable, sympy.together(_ilt(expr))))
+
+        return _ilt(e)
+
+    assert M.is_square
+    N = M.shape[0]
+    s = sympy.Dummy('s')
+
+    Ms = (s*sympy.eye(N) - M)
+    Mres = Ms.adjugate() / Ms.det()
+
+    def expMij(i, j):
+        """Partial fraction expansion then inverse Laplace transform"""
+        Mresij_pfe = sympy.apart(Mres[i, j], s, full=True)
+        return ilt(Mresij_pfe, s, t)
+
+    return sympy.Matrix(N, N, expMij)
+
+
diff --git a/odetoolbox/system_of_shapes.py b/odetoolbox/system_of_shapes.py
index 164b77ab..a41cda43 100644
--- a/odetoolbox/system_of_shapes.py
+++ b/odetoolbox/system_of_shapes.py
@@ -32,7 +32,7 @@
 from .config import Config
 from .shapes import Shape
 from .singularity_detection import SingularityDetection, SingularityDetectionException
-from .sympy_helpers import _custom_simplify_expr, _is_zero
+from .sympy_helpers import _custom_simplify_expr, _is_zero, expMt
 
 
 class GetBlockDiagonalException(Exception):
@@ -46,6 +46,7 @@ def get_block_diagonal_blocks(A):
     assert A.shape[0] == A.shape[1], "matrix A should be square"
 
     A_connectivity_undirected = (A != 0) | (A.T != 0)    # make symmetric (undirected) connectivity graph from the system matrix
+    A_connectivity_undirected = np.triu(A_connectivity_undirected)
 
     graph_components = scipy.sparse.csgraph.connected_components(A_connectivity_undirected)[1]
 
@@ -88,7 +89,6 @@ def __init__(self, x: sympy.Matrix, A: sympy.Matrix, b: sympy.Matrix, c: sympy.M
         :param b: Vector containing inhomogeneous part (constant term).
         :param c: Vector containing nonlinear part.
         """
-        logging.debug("Initializing system of shapes with x = " + str(x) + ", A = " + str(A) + ", b = " + str(b) + ", c = " + str(c))
         assert x.shape[0] == A.shape[0] == A.shape[1] == b.shape[0] == c.shape[0]
         self.x_ = x
         self.A_ = A
@@ -214,7 +214,13 @@ def _generate_propagator_matrix(self, A):
         # optimized: be explicit about block diagonal elements; much faster!
         try:
             blocks = get_block_diagonal_blocks(np.array(A))
-            propagators = [sympy.simplify(sympy.exp(sympy.Matrix(block) * sympy.Symbol(Config().output_timestep_symbol))) for block in blocks]
+
+            if Config().use_alternative_expM:
+                expM = expMt
+            else:
+                expM = sympy.exp
+
+            propagators = [sympy.simplify(expM(sympy.Matrix(block) * sympy.Symbol(Config().output_timestep_symbol, real=True))) for block in blocks]
             P = sympy.Matrix(scipy.linalg.block_diag(*propagators))
         except GetBlockDiagonalException:
             # naive: calculate propagators in one step
@@ -244,12 +250,12 @@ def generate_propagator_solver(self, disable_singularity_detection: bool = False
                 if conditions:
                     # if there is one or more condition under which the solution goes to infinity...
 
-                    logging.warning("Under certain conditions, the propagator matrix is singular (contains infinities).")
-                    logging.warning("List of all conditions that result in a division by zero:")
+                    logging.getLogger(__name__).warning("Under certain conditions, the propagator matrix is singular (contains infinities).")
+                    logging.getLogger(__name__).warning("List of all conditions that result in a division by zero:")
                     for cond_set in conditions:
-                        logging.warning("\t" + r" ∧ ".join([str(eq.lhs) + " = " + str(eq.rhs) for eq in cond_set]))
+                        logging.getLogger(__name__).warning("\t" + r" ∧ ".join([str(eq.lhs) + " = " + str(eq.rhs) for eq in cond_set]))
             except SingularityDetectionException:
-                logging.warning("Could not check the propagator matrix for singularities.")
+                logging.getLogger(__name__).warning("Could not check the propagator matrix for singularities.")
 
         #
         #   generate symbols for each nonzero entry of the propagator matrix
@@ -305,7 +311,6 @@ def generate_propagator_solver(self, disable_singularity_detection: bool = False
             if not _is_zero(self.b_[row]):
                 # only simplify in case an inhomogeneous term is present
                 update_expr[str(self.x_[row])] = _custom_simplify_expr(update_expr[str(self.x_[row])])
-            logging.debug("update_expr[" + str(self.x_[row]) + "] = " + str(update_expr[str(self.x_[row])]))
 
         all_state_symbols = [str(sym) for sym in self.x_]
         initial_values = {sym: str(self.get_initial_value(sym)) for sym in all_state_symbols}