diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0b14fbc..64ac6f8 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -110,13 +110,11 @@ include_directories(
   src/sys
   )
 
-link_directories(/usr/lib/chaiscript)
-
 # Add the code subdirectory
 add_subdirectory(src)
 add_subdirectory(testing)
 
 add_executable(tasmet src/main.cpp src/gui/tasmet_resources.qrc)
 add_executable(tasmet_solvemodel src/tasmet_solvemodel.cpp)
-target_link_libraries(tasmet tasmet_gui tasmet_src messages Qt5::Widgets chaiscript_stdlib-5.7.0 openblas)
-target_link_libraries(tasmet_solvemodel tasmet_src messages openblas chaiscript_stdlib-5.7.0)
+target_link_libraries(tasmet tasmet_gui tasmet_src messages Qt5::Widgets muparser openblas)
+target_link_libraries(tasmet_solvemodel tasmet_src messages openblas muparser)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index f0b4a7b..0906a5c 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -27,7 +27,6 @@ include_directories(
   )
 
 add_library(tasmet_src
-  chaiscript.cpp
 
   tasmet_tracer.cpp
   tasmet_exception.cpp
diff --git a/src/chaiscript.cpp b/src/chaiscript.cpp
deleted file mode 100644
index 4ab5873..0000000
--- a/src/chaiscript.cpp
+++ /dev/null
@@ -1,28 +0,0 @@
-// chaiscript.cpp
-//
-// last-edit-by: J.A. de Jong 
-// 
-// Description:
-//
-//////////////////////////////////////////////////////////////////////
-
-#include "chaiscript.h"
-#include <chaiscript/chaiscript_stdlib.hpp>
-#include <chaiscript/chaiscript.hpp>
-#include "chaiscript/math.hpp"
-#include "tasmet_types.h"
-
-using chaiscript::ChaiScript;
-
-std::unique_ptr<chaiscript::ChaiScript> getChaiScriptInstance() {
-
-    auto mathlib = chaiscript::extras::math::bootstrap();
-    
-    std::unique_ptr<ChaiScript> chai(new ChaiScript());
-    // std::unique_ptr<ChaiScript> chai(new ChaiScript(chaiscript::Std_Lib::library()));
-    
-    chai->add(mathlib);
-    chai->add_global(chaiscript::var(number_pi),"pi");
-    return chai;
-}
-//////////////////////////////////////////////////////////////////////
diff --git a/src/chaiscript.h b/src/chaiscript.h
deleted file mode 100644
index 80877c7..0000000
--- a/src/chaiscript.h
+++ /dev/null
@@ -1,20 +0,0 @@
-// chaiscript.h
-//
-// Author: J.A. de Jong 
-//
-// Description:
-// Slow compiling chaiscript therefore we provide this wrapper
-//////////////////////////////////////////////////////////////////////
-#pragma once
-#ifndef CHAISCRIPT_H
-#define CHAISCRIPT_H
-namespace chaiscript {
-    class ChaiScript;
-}
-#include <memory>
-
-std::unique_ptr<chaiscript::ChaiScript> getChaiScriptInstance();
-
-#endif // CHAISCRIPT_H
-//////////////////////////////////////////////////////////////////////
-
diff --git a/src/chaiscript/math.hpp b/src/chaiscript/math.hpp
deleted file mode 100644
index 9d5b2a3..0000000
--- a/src/chaiscript/math.hpp
+++ /dev/null
@@ -1,809 +0,0 @@
-#include <cmath>
-#include <memory>
-
-#include <chaiscript/chaiscript.hpp>
-
-namespace chaiscript {
-  namespace extras {
-    namespace math {
-      // TRIG FUNCTIONS
-      template<typename Ret, typename Param>
-      ModulePtr cos(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::cos)), "cos");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr sin(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::sin)), "sin");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr tan(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::tan)), "tan");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr acos(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::acos)), "acos");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr asin(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::asin)), "asin");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr atan(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::atan)), "atan");
-        return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr atan2(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::atan2)), "atan2");
-        return m;
-      }
-
-      // HYPERBOLIC FUNCTIONS
-      template<typename Ret, typename Param>
-      ModulePtr cosh(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::cosh)), "cosh");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr sinh(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::sinh)), "sinh");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr tanh(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::tanh)), "tanh");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr acosh(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::acosh)), "acosh");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr asinh(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::asinh)), "asinh");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr atanh(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::atanh)), "atanh");
-        return m;
-      }
-
-      // EXPONENTIAL AND LOGARITHMIC FUNCTIONS
-      template<typename Ret, typename Param>
-      ModulePtr exp(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::exp)), "exp");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr frexp(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::frexp)), "frexp");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr ldexp(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::ldexp)), "ldexp");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr log(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::log)), "log");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr log10(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::log10)), "log10");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr modf(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::modf)), "modf");
-    	  return m;
-      }
-      template<typename Ret, typename Param>
-      ModulePtr exp2(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::exp2)), "exp2");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr expm1(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::expm1)), "expm1");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr ilogb(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::ilogb)), "ilogb");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr log1p(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::log1p)), "log1p");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr log2(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::log2)), "log2");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr logb(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::logb)), "logb");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr scalbn(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::scalbn)), "scalbn");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr scalbln(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::scalbln)), "scalbln");
-    	  return m;
-      }
-
-      // POWER FUNCTIONS
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr pow(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::pow)), "pow");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr sqrt(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::sqrt)), "sqrt");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr cbrt(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::cbrt)), "cbrt");
-        return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr hypot(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::hypot)), "hypot");
-        return m;
-      }
-
-      // ERROR AND GAMMA FUNCTIONS
-      template<typename Ret, typename Param>
-      ModulePtr erf(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::erf)), "erf");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr erfc(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::erfc)), "erfc");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr tgamma(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::tgamma)), "tgamma");
-        return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr lgamma(ModulePtr m = std::make_shared<Module>())
-      {
-        m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::lgamma)), "lgamma");
-        return m;
-      }
-
-      // ROUNDING AND REMAINDER FUNCTIONS
-      template<typename Ret, typename Param>
-      ModulePtr ceil(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::ceil)), "ceil");
-		  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr floor(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::floor)), "floor");
-		  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr fmod(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::fmod)), "fmod");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr trunc(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::trunc)), "trunc");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr round(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::round)), "round");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr lround(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::lround)), "lround");
-    	  return m;
-      }
-
-      // long long ints do not work
-      template<typename Ret, typename Param>
-      ModulePtr llround(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::llround)), "llround");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr rint(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::rint)), "rint");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr lrint(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::lrint)), "lrint");
-    	  return m;
-      }
-
-      // long long ints do not work
-      template<typename Ret, typename Param>
-      ModulePtr llrint(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::llrint)), "llrint");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr nearbyint(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::nearbyint)), "nearbyint");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr remainder(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::remainder)), "remainder");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2, typename Param3>
-      ModulePtr remquo(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2, Param3)>(&std::remquo)), "remquo");
-    	  return m;
-      }
-
-      // FLOATING-POINT MANIPULATION FUNCTIONS
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr copysign(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::copysign)), "copysign");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr nan(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::nan)), "nan");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr nextafter(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::nextafter)), "nextafter");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr nexttoward(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::nexttoward)), "nexttoward");
-    	  return m;
-      }
-
-      // MINIMUM, MAXIMUM, DIFFERENCE FUNCTIONS
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr fdim(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::fdim)), "fdim");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr fmax(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::fmax)), "fmax");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr fmin(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::fmin)), "fmin");
-    	  return m;
-      }
-
-      // OTHER FUNCTIONS
-      template<typename Ret, typename Param>
-      ModulePtr fabs(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::fabs)), "fabs");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr abs(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::abs)), "abs");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2, typename Param3>
-      ModulePtr fma(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2, Param3)>(&std::fma)), "fma");
-    	  return m;
-      }
-
-      // CLASSIFICATION FUNCTIONS
-      template<typename Ret, typename Param>
-      ModulePtr fpclassify(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::fpclassify)), "fpclassify");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr isfinite(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::isfinite)), "isfinite");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr isinf(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::isinf)), "isinf");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr isnan(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::isnan)), "isnan");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr isnormal(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::isnormal)), "isnormal");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param>
-      ModulePtr signbit(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::signbit)), "signbit");
-    	  return m;
-      }
-
-
-      // COMPARISON FUNCTIONS
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr isgreater(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::isgreater)), "isgreater");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr isgreaterequal(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::isgreaterequal)), "isgreaterequal");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr isless(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::isless)), "isless");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr islessequal(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::islessequal)), "islessequal");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr islessgreater(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::islessgreater)), "islessgreater");
-    	  return m;
-      }
-
-      template<typename Ret, typename Param1, typename Param2>
-      ModulePtr isunordered(ModulePtr m = std::make_shared<Module>())
-      {
-    	  m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::isunordered)), "isunordered");
-    	  return m;
-      }
-
-      ModulePtr bootstrap(ModulePtr m = std::make_shared<Module>())
-      {
-    	// TRIG FUNCTIONS
-        cos<double, double>(m);
-        cos<float, float>(m);
-        cos<long double, long double>(m);
-
-        sin<double, double>(m);
-        sin<float, float>(m);
-        sin<long double, long double>(m);
-
-        tan<double, double>(m);
-        tan<float, float>(m);
-        tan<long double, long double>(m);
-
-        acos<double, double>(m);
-        acos<float, float>(m);
-        acos<long double, long double>(m);
-
-        asin<double, double>(m);
-        asin<float, float>(m);
-        asin<long double, long double>(m);
-
-        atan<double, double>(m);
-        atan<float, float>(m);
-        atan<long double, long double>(m);
-
-        atan2<double, double, double>(m);
-        atan2<float, float, float>(m);
-        atan2<long double, long double, long double>(m);
-
-        // HYPERBOLIC FUNCTIONS
-        cosh<double, double>(m);
-        cosh<float, float>(m);
-        cosh<long double, long double>(m);
-
-        sinh<double, double>(m);
-        sinh<float, float>(m);
-        sinh<long double, long double>(m);
-
-        tanh<double, double>(m);
-        tanh<float, float>(m);
-        tanh<long double, long double>(m);
-
-        acosh<double, double>(m);
-        acosh<float, float>(m);
-        acosh<long double, long double>(m);
-
-        asinh<double, double>(m);
-        asinh<float, float>(m);
-        asinh<long double, long double>(m);
-
-        atanh<double, double>(m);
-        atanh<float, float>(m);
-        atanh<long double, long double>(m);
-
-        // EXPONENTIAL AND LOGARITHMIC FUNCTIONS
-        exp<double, double>(m);
-        exp<float, float>(m);
-        exp<long double, long double>(m);
-
-        frexp<double, double, int *>(m);
-        frexp<float, float, int *>(m);
-        frexp<long double, long double, int *>(m);
-
-        ldexp<double, double, int>(m);
-        ldexp<float, float, int>(m);
-        ldexp<long double, long double, int>(m);
-
-        log<double, double>(m);
-        log<float, float>(m);
-        log<long double, long double>(m);
-
-        log10<double, double>(m);
-        log10<float, float>(m);
-        log10<long double, long double>(m);
-
-        modf<double, double, double *>(m);
-        modf<float, float, float *>(m);
-        modf<long double, long double, long double *>(m);
-
-        exp2<double, double>(m);
-        exp2<float, float>(m);
-        exp2<long double, long double>(m);
-
-        expm1<double, double>(m);
-        expm1<float, float>(m);
-        expm1<long double, long double>(m);
-
-        ilogb<int, double>(m);
-        ilogb<int, float>(m);
-        ilogb<int, long double>(m);
-
-        log1p<double, double>(m);
-        log1p<float, float>(m);
-        log1p<long double, long double>(m);
-
-        log2<double, double>(m);
-        log2<float, float>(m);
-        log2<long double, long double>(m);
-
-        logb<double, double>(m);
-        logb<float, float>(m);
-        logb<long double, long double>(m);
-
-        scalbn<double, double, int>(m);
-        scalbn<float, float, int>(m);
-        scalbn<long double, long double, int>(m);
-
-        scalbln<double, double, long int>(m);
-        scalbln<float, float, long int>(m);
-        scalbln<long double, long double, long int>(m);
-
-        // POWER FUNCTIONS
-        pow<double, double, double>(m);
-        pow<float, float, float>(m);
-        pow<long double, long double, long double>(m);
-
-        sqrt<double, double>(m);
-        sqrt<float, float>(m);
-        sqrt<long double, long double>(m);
-
-        cbrt<double, double>(m);
-        cbrt<float, float>(m);
-        cbrt<long double, long double>(m);
-
-        hypot<double, double, double>(m);
-        hypot<float, float, float>(m);
-        hypot<long double, long double, long double>(m);
-
-        // ERROR AND GAMMA FUNCTIONS
-        erf<double, double>(m);
-        erf<float, float>(m);
-        erf<long double, long double>(m);
-
-        erfc<double, double>(m);
-        erfc<float, float>(m);
-        erfc<long double, long double>(m);
-
-        tgamma<double, double>(m);
-        tgamma<float, float>(m);
-        tgamma<long double, long double>(m);
-
-        lgamma<double, double>(m);
-        lgamma<float, float>(m);
-        lgamma<long double, long double>(m);
-
-        // ROUNDING AND REMAINDER FUNCTIONS
-        ceil<double, double>(m);
-        ceil<float, float>(m);
-        ceil<long double, long double>(m);
-
-        floor<double, double>(m);
-        floor<float, float>(m);
-        floor<long double, long double>(m);
-
-        fmod<double, double, double>(m);
-        fmod<float, float, float>(m);
-        fmod<long double, long double, long double>(m);
-
-        trunc<double, double>(m);
-        trunc<float, float>(m);
-        trunc<long double, long double>(m);
-
-        round<double, double>(m);
-        round<float, float>(m);
-        round<long double, long double>(m);
-
-        lround<long int, double>(m);
-        lround<long int, float>(m);
-        lround<long int, long double>(m);
-
-        // long long ints do not work
-        llround<long long int, double>(m);
-        llround<long long int, float>(m);
-        llround<long long int, long double>(m);
-
-        rint<double, double>(m);
-        rint<float, float>(m);
-        rint<long double, long double>(m);
-
-        lrint<long int, double>(m);
-        lrint<long int, float>(m);
-        lrint<long int, long double>(m);
-
-        // long long ints do not work
-        llrint<long long int, double>(m);
-        llrint<long long int, float>(m);
-        llrint<long long int, long double>(m);
-
-        nearbyint<double, double>(m);
-        nearbyint<float, float>(m);
-        nearbyint<long double, long double>(m);
-
-        remainder<double, double, double>(m);
-        remainder<float, float, float>(m);
-        remainder<long double, long double, long double>(m);
-
-        remquo<double, double, double, int *>(m);
-        remquo<float, float, float, int *>(m);
-        remquo<long double, long double, long double, int *>(m);
-
-        // FLOATING-POINT MANIPULATION FUNCTIONS
-        copysign<double, double, double>(m);
-        copysign<float, float, float>(m);
-        copysign<long double, long double, long double>(m);
-
-        nan<double, const char*>(m);
-
-        nextafter<double, double, double>(m);
-        nextafter<float, float, float>(m);
-        nextafter<long double, long double, long double>(m);
-
-        nexttoward<double, double, long double>(m);
-        nexttoward<float, float, long double>(m);
-        nexttoward<long double, long double, long double>(m);
-
-        // MINIMUM, MAXIMUM, DIFFERENCE FUNCTIONS
-        fdim<double, double, double>(m);
-        fdim<float, float, float>(m);
-        fdim<long double, long double, long double>(m);
-
-        fmax<double, double, double>(m);
-        fmax<float, float, float>(m);
-        fmax<long double, long double, long double>(m);
-
-        fmin<double, double, double>(m);
-        fmin<float, float, float>(m);
-        fmin<long double, long double, long double>(m);
-
-        // OTHER FUNCTIONS
-        fabs<double, double>(m);
-        fabs<float, float>(m);
-        fabs<long double, long double>(m);
-
-        abs<double, double>(m);
-        abs<float, float>(m);
-        abs<long double, long double>(m);
-
-        fma<double, double, double, double>(m);
-        fma<float, float, float, float>(m);
-        fma<long double, long double, long double, long double>(m);
-
-        // CLASSIFICATION FUNCTIONS
-        fpclassify<int, float>(m);
-        fpclassify<int, double>(m);
-        fpclassify<int, long double>(m);
-
-        isfinite<bool, float>(m);
-        isfinite<bool, double>(m);
-        isfinite<bool, long double>(m);
-
-        isinf<bool, float>(m);
-        isinf<bool, double>(m);
-        isinf<bool, long double>(m);
-
-        isnan<bool, float>(m);
-        isnan<bool, double>(m);
-        isnan<bool, long double>(m);
-
-        isnormal<bool, float>(m);
-        isnormal<bool, double>(m);
-        isnormal<bool, long double>(m);
-
-        signbit<bool, float>(m);
-        signbit<bool, double>(m);
-        signbit<bool, long double>(m);
-
-        // COMPARISON FUNCTIONS
-        isgreater<bool, double, double>(m);
-        isgreater<bool, float, float>(m);
-        isgreater<bool, long double, long double>(m);
-
-        isgreaterequal<bool, double, double>(m);
-        isgreaterequal<bool, float, float>(m);
-        isgreaterequal<bool, long double, long double>(m);
-
-        isless<bool, double, double>(m);
-        isless<bool, float, float>(m);
-        isless<bool, long double, long double>(m);
-
-        islessequal<bool, double, double>(m);
-        islessequal<bool, float, float>(m);
-        islessequal<bool, long double, long double>(m);
-
-        islessgreater<bool, double, double>(m);
-        islessgreater<bool, float, float>(m);
-        islessgreater<bool, long double, long double>(m);
-
-        isunordered<bool, double, double>(m);
-        isunordered<bool, float, float>(m);
-        isunordered<bool, long double, long double>(m);
-
-        return m;
-      }
-    }
-  }
-}
diff --git a/src/duct/duct.cpp b/src/duct/duct.cpp
index fc39f21..da0f6f2 100644
--- a/src/duct/duct.cpp
+++ b/src/duct/duct.cpp
@@ -346,7 +346,7 @@ void Duct::exportHDF5(const hid_t group_id) const {
 
     PosData S;
     S.name = "Cross-sectional area";
-    S.unit = "m^2";
+    S.unit = "m$^2$";
     S.symbol = "S";
     S.x = this->S;
     S.exportHDF5(group_id);
diff --git a/src/tasmet_constants.h b/src/tasmet_constants.h
index d136e52..1c71748 100644
--- a/src/tasmet_constants.h
+++ b/src/tasmet_constants.h
@@ -14,25 +14,25 @@
 // Unfortunately to let the code compile with Swig v 2.0, strongly
 // typed enums are not supported. Therefore this is a normal
 // enumerated type and not an enum class.
-DECLARE_ENUM(Varnr,
-             none,                         // None
-             rho,                        // Density
-             m,                          // Mass flow (rho*U)
-             T,                          // Temperature
-             p,                          // Pressure
-             Ts,                         // Temperature of the solid
-             Tw,                         // Temperature of the solid wall
-             mH,                         // Enthalpy flow (Watts)
-             U,                          // Volume flow (m^3/s)
-             u,                          // Velocity (U/Sf)
-             mu,                         // Momentum flux
-             Q,                          // Heat flow
-             Qs,                 // Solid heat Flow
-             F,                 // A mechanical domain force [N]
-             x,                // A mechanical displacement [m]
-             Z,                 // A mechanical impedance [N/m]
-             mEkin	       // Kinetic energy flow (Watts)
-             );
+// DECLARE_ENUM(Varnr,
+//              none,                         // None
+//              rho,                        // Density
+//              m,                          // Mass flow (rho*U)
+//              T,                          // Temperature
+//              p,                          // Pressure
+//              Ts,                         // Temperature of the solid
+//              Tw,                         // Temperature of the solid wall
+//              mH,                         // Enthalpy flow (Watts)
+//              U,                          // Volume flow (m^3/s)
+//              u,                          // Velocity (U/Sf)
+//              mu,                         // Momentum flux
+//              Q,                          // Heat flow
+//              Qs,                 // Solid heat Flow
+//              F,                 // A mechanical domain force [N]
+//              x,                // A mechanical displacement [m]
+//              Z,                 // A mechanical impedance [N/m]
+//              mEkin	       // Kinetic energy flow (Watts)
+//              );
 
 DECLARE_ENUM(Pos,posleft,posright);
 
diff --git a/src/tasmet_evalscript.cpp b/src/tasmet_evalscript.cpp
index 494f1d2..ff2aae8 100644
--- a/src/tasmet_evalscript.cpp
+++ b/src/tasmet_evalscript.cpp
@@ -5,85 +5,85 @@
 // Description:
 //
 //////////////////////////////////////////////////////////////////////
-
+// #define TRACERPLUS 16
 #include "tasmet_evalscript.h"
-#include <chaiscript/chaiscript.hpp>
-#include "chaiscript.h"
 #include "tasmet_exception.h"
 #include "tasmet_tracer.h"
 #include <sstream>
-using chaiscript::ChaiScript;
-using std::stringstream;
+#include <muParser.h>
 
-template<typename T>
-inline T wrap_eval(ChaiScript* chai,const string& script) {
-   try {
-       return chai->eval<T>(script);
-    }
-    catch(chaiscript::exception::eval_error &e) {
-        TRACE(15,script);
-        throw TaSMETError(e.what());
-    }
-   catch(std::bad_cast &e) {
-        TRACE(15,script);
-        throw TaSMETError("Cannot get return value from script");
-    }
-
-}
-template<>
-inline void wrap_eval<void>(ChaiScript* chai,const string& script) {
-   try {
-       chai->eval(script);
-    }
-   catch(std::runtime_error &e) {
-        TRACE(15,script);
-        throw TaSMETError(e.what());
-    }
-   catch(std::bad_cast &e) {
-        TRACE(15,script);
-        throw TaSMETError("Cannot get return value from script");
-    }
-
-}
+const int vector_prealloc = 30;
 
 EvaluateFun::EvaluateFun(const string& fun_return,
                          const string& err_msg,
-                         const string& vars):
+                         const string& var):
     _err_msg(err_msg),
     _fun_return(fun_return)
 {
-    TRACE(15,"EvaluateFun::EvaluateFun()");
-    _chai = getChaiScriptInstance();
-    if(!_chai) throw TaSMETBadAlloc();
 
-    string script = string("def myfun(") + vars + ") {\n";
-    script += "return " + fun_return + "; }\n";
-    
-    wrap_eval<void>(_chai.get(),script);
+    TRACE(15,"EvaluateFun::EvaluateFun()");
+
+    _parser = new mu::Parser();
+    if(!_parser) throw TaSMETBadAlloc();
+
+    /**
+     * It is imperative to put a reserve operation here, which
+     * preallocates enough space in the vector. Otherwise a
+     * reallocation may invalidate the pointers which the parser
+     * points to.
+     * 
+     */
+    _values.reserve(vector_prealloc);
+
+    _values.push_back(0);       // This is the (x), or (t)
+
+    // Add pi, as it unfortunately does not exist
+    _values.push_back(number_pi);
+    try {
+        _parser->DefineVar(var,&_values[0]);
+        _parser->DefineVar("pi",&_values[1]);
+        _parser->SetExpr(fun_return);
+    }
+    catch(mu::Parser::exception_type& e) {
+        TRACE(15,e.GetMsg());
+        throw TaSMETError(_err_msg);
+    }
     
 }
 void EvaluateFun::addGlobalDef(const string& name,const d value) {
     TRACE(15,"EvaluateFun::addGlobalDef()");
-
-    _chai->add_global(chaiscript::var(value),name);
+    if(_values.size() == vector_prealloc -1)
+        throw TaSMETError("Too many definitions");
+    
+    _values.push_back(value);
+    try {    
+        _parser->DefineVar(name,&_values[_values.size()-1]);
+    }
+    catch(mu::Parser::exception_type& e) {
+        TRACE(15,e.GetMsg());
+        throw TaSMETError(_err_msg);
+    }
 }
 vd EvaluateFun::operator()(const vd& x) {
     TRACE(15,"EvaluateFun::operator(vd)");
     vd y(x.size());
 
     for(us i=0;i<x.size();i++) {
-
-        stringstream value;
-        value << "myfun(" << std::scientific << x(i) << ");";
-
-        y(i) = wrap_eval<d>(_chai.get(),value.str());
-        
+        _values[0] = x[i];
+        VARTRACE(15,_values[0]);
+        try {
+            y(i) = _parser->Eval();
+        }
+        catch(mu::Parser::exception_type& e) {
+            TRACE(15,e.GetMsg());
+            throw TaSMETError(_err_msg);
+        }
     }
     return y;
 }
 
 EvaluateFun::~EvaluateFun() {
-
+    delete _parser;
 }
 
 
diff --git a/src/tasmet_evalscript.h b/src/tasmet_evalscript.h
index da26eaf..b668969 100644
--- a/src/tasmet_evalscript.h
+++ b/src/tasmet_evalscript.h
@@ -13,24 +13,50 @@
 #include "tasmet_types.h"
 #include <memory>
 
-namespace chaiscript {
-    class ChaiScript;
+namespace mu {
+    class Parser;
 }
-
 // Helper class to evaluate simple 1D math function evaluations in Python
 class EvaluateFun {
-    std::unique_ptr<chaiscript::ChaiScript> _chai;
     string _err_msg;
     string _fun_return;
+    mu::Parser* _parser = nullptr;
+    std::vector<d> _values;
 public:
-    EvaluateFun(const string& fun_return,const string& err_msg = "Script error",const string& vars = "x");
+    /** 
+     * Construct a new function evaluator
+     *
+     * @param fun_return what the function should return
+     * (e.g. sin(x)).  
+     *
+     * @param err_msg The exception message that
+     * should be thrown in case something goes wrong.
+     *
+     * @param var The variable on which the function depends (e.g. x,
+     * or t)
+     */
+    EvaluateFun(const string& fun_return,
+                const string& err_msg = "Script error",
+                const string& var = "x");
     
-    // Add a global definition to the namespace
+    /** 
+     * Add a global parameter to the namespace.
+     *
+     * @param name 
+     * @param value 
+     */
     void addGlobalDef(const string& name,
                       const d value);
     
-    // Evaluate a single function at multiple points and return
-    // the result for each point
+    /** 
+     * Evaluate the  function at multiple points and return
+     * the result for each point
+     * 
+     * @param points: the x-values, or t-values where to evaluate the
+     * function.
+     *
+     * @return the results of the function evalated at the points
+     */
     vd operator()(const vd& points);
 
     // Used to cleanup the python namespace