diff --git a/.gitignore b/.gitignore
index 346a53f..0e398a2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -21,3 +21,4 @@ doc/usg.pdf
 brent_test
 tasmet
 tasmet_automoc.dir
+tasmet_solvemodel
diff --git a/CMakeLists.txt b/CMakeLists.txt
index ef0a6c1..1b559af 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -126,11 +126,14 @@ include_directories(
   src/material
   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 PythonQt Qt5::Widgets openblas tcmalloc)
-target_link_libraries(tasmet_solvemodel tasmet_src messages PythonQt openblas)
+target_link_libraries(tasmet tasmet_gui tasmet_src messages PythonQt Qt5::Widgets chaiscript_stdlib-5.7.0 openblas)
+target_link_libraries(tasmet_solvemodel tasmet_src messages PythonQt openblas chaiscript_stdlib-5.7.0)
diff --git a/Doxyfile b/Doxyfile
index 515149c..dad1477 100644
--- a/Doxyfile
+++ b/Doxyfile
@@ -187,7 +187,7 @@ SHORT_NAMES            = NO
 # description.)
 # The default value is: NO.
 
-JAVADOC_AUTOBRIEF      = NO
+JAVADOC_AUTOBRIEF      = YES
 
 # If the QT_AUTOBRIEF tag is set to YES then doxygen will interpret the first
 # line (until the first dot) of a Qt-style comment as the brief description. If
diff --git a/doc/usg.lyx b/doc/usg.lyx
index cebbed8..3d6d3bd 100644
--- a/doc/usg.lyx
+++ b/doc/usg.lyx
@@ -439,19 +439,6 @@ Solver
  object.
  Often, the user creates a Python scripts which has approximately the following
  shape:
-\begin_inset ERT
-status open
-
-\begin_layout Plain Layout
-
-
-\backslash
-lstinputlisting{tutorial/program_overview.py}
-\end_layout
-
-\end_inset
-
-
 \end_layout
 
 \begin_layout Standard
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index e5f709f..632f63a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -27,10 +27,12 @@ include_directories(
   )
 
 add_library(tasmet_src
+  chaiscript.cpp
+
   tasmet_tracer.cpp
   tasmet_exception.cpp
   tasmet_assert.cpp
-  tasmet_pyeval.cpp
+  tasmet_evalscript.cpp
 
   protobuf/message_tools.cpp
 
@@ -72,4 +74,4 @@ add_library(tasmet_src
 
 add_subdirectory(gui) 
 add_subdirectory(protobuf)
-target_link_libraries(tasmet_src Qt5::Widgets)
+target_link_libraries(tasmet_src Qt5::Widgets superlu dl pthread)
diff --git a/src/chaiscript.cpp b/src/chaiscript.cpp
new file mode 100644
index 0000000..4ab5873
--- /dev/null
+++ b/src/chaiscript.cpp
@@ -0,0 +1,28 @@
+// 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
new file mode 100644
index 0000000..80877c7
--- /dev/null
+++ b/src/chaiscript.h
@@ -0,0 +1,20 @@
+// 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
new file mode 100644
index 0000000..9d5b2a3
--- /dev/null
+++ b/src/chaiscript/math.hpp
@@ -0,0 +1,809 @@
+#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 2fbae1e..a978296 100644
--- a/src/duct/duct.cpp
+++ b/src/duct/duct.cpp
@@ -8,24 +8,27 @@
 #include "duct.h"
 #include "tasystem.h"
 #include "tasmet_assert.h"
-#include "tasmet_pyeval.h"
+#include "tasmet_evalscript.h"
 
-Duct::Duct(const us id,const pb::Duct& duct):
-    Segment(id,duct.name()),
-    Geom(duct)
+
+Duct::Duct(const us id,const pb::Duct& ductpb):
+    Segment(id,ductpb.name()),
+    Geom(ductpb),
+    _ductpb(ductpb)
 {
     TRACE(15,"Duct::Duct()");
+
     const char* invTsfun = "Invalid solid-temperature prescribing function";
-    EvaluateFun Tsfun(duct.stempfunc(),invTsfun);
-    Tsfun.addGlobalDef("L",duct.length());
+    EvaluateFun Tsfun(ductpb.stempfunc(),invTsfun);
+    Tsfun.addGlobalDef("L",ductpb.length());
     _Tsprescribed = Tsfun(x);
+
     if(min(_Tsprescribed) < constants::min_T0 ||
        max(_Tsprescribed) > constants::max_T0) {
         throw TaSMETError(invTsfun);
     }
 
-
-    switch (duct.htmodel()) {
+    switch (ductpb.htmodel()) {
     case pb::Isentropic: {
         break;
     }
@@ -34,25 +37,31 @@ Duct::Duct(const us id,const pb::Duct& duct):
     break;
     }
 
-
 }
 Duct::Duct(const Duct& other):
     Segment(other),
     Geom(other),
+    _ductpb(other._ductpb),
     _Tsprescribed(other._Tsprescribed)
 {
     // Do something with the equations here
+    TRACE(15,"Duct::~Duct");
+
+    
+
 }
 Duct* Duct::copy() const {
     return new Duct(*this);
 }
 Duct::~Duct() {
+    TRACE(15,"Duct::~Duct");
     // for(Equation* eq: _eqs){
     //     delete eq;
     // }
 }
 void Duct::residual(const TaSystem& sys,arma::subview_col<d> && residual) const {
 
+    TRACE(15,"Duct::residual()");
     const arma::subview_col<d> sol = sys.getSolution(_id);
 
     VARTRACE(15,sol(0));
@@ -61,16 +70,22 @@ void Duct::residual(const TaSystem& sys,arma::subview_col<d> && residual) const
 }
 vd Duct::initialSolution(const TaSystem& sys) const {
     
-    vd initsol(getNDofs(sys));
+    TRACE(15,"Duct::initialSolution()");
+
+    vd initsol(getNDofs(sys));
+    VARTRACE(15,initsol.size());
+
+    us vars_per_gp = 4;
+    vars_per_gp+= (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
 
-    us vars_per_gp = 5;
     us Ns = sys.Ns();
 
     const Gas& gas = sys.gas();
     
     for(us i=0;i<ngp();i++){
-
+        VARTRACE(15,i);
         // Initial density
+
         initsol.subvec((i*vars_per_gp+0)*Ns,(i*vars_per_gp+1)*Ns-1) =
             gas.rho0();
         
@@ -79,17 +94,18 @@ vd Duct::initialSolution(const TaSystem& sys) const {
             0;
 
         // Initial Temperature
-        initsol.subvec((i*vars_per_gp+2)*Ns,(i*vars_per_gp+2)*Ns-1) =
+        initsol.subvec((i*vars_per_gp+2)*Ns,(i*vars_per_gp+3)*Ns-1) =
             _Tsprescribed(i);
 
         // Initial pressure
         initsol.subvec((i*vars_per_gp+3)*Ns,(i*vars_per_gp+4)*Ns-1) =
             gas.p0();
         
-        // Initial solid temperature
-        initsol.subvec((i*vars_per_gp+4)*Ns,(i*vars_per_gp+5)*Ns-1) =
-            _Tsprescribed(i);
-
+        // Initial solid temperature, if not prescribed
+        if(_ductpb.stempmodel() != pb::Prescribed) {
+            initsol.subvec((i*vars_per_gp+4)*Ns,(i*vars_per_gp+5)*Ns-1) =
+                _Tsprescribed(i);
+        }
     }
 
     return initsol;
@@ -97,18 +113,23 @@ vd Duct::initialSolution(const TaSystem& sys) const {
 }
 
 us Duct::getNEqs(const TaSystem& sys) const {
-
+    TRACE(15,"Duct::getNEqs()");    
     us Ns = sys.Ns();
-    us number_eqs = _eqs.size();
+    us number_eqs = 4;
+
+    // When we have to solve a solid heat balance
+    number_eqs+= (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
+
+    return Ns*number_eqs*(ngp()-1);
 
-    return Ns*number_eqs*ngp();
 }
 us Duct::getNDofs(const TaSystem& sys) const {
-    
+    TRACE(15,"Duct::getNDofs()");        
     us Ns = sys.Ns();  
 
     // rho,u,T,p,Ts
-    us nvars_per_gp = 5;
+    us nvars_per_gp = 4;
+    nvars_per_gp += (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
 
     return Ns*nvars_per_gp*ngp();
 }
diff --git a/src/duct/duct.h b/src/duct/duct.h
index c036087..4bad365 100644
--- a/src/duct/duct.h
+++ b/src/duct/duct.h
@@ -22,8 +22,11 @@ class Duct : public Segment, public Geom {
     // Drag* _drag = nullptr;
     // Heat* _heat = nullptr;
     std::vector<Equation*> _eqs;
+    pb::Duct _ductpb;
     vd _Tsprescribed;
 
+
+
 protected:
     Duct(const Duct&);
 public:
@@ -38,6 +41,8 @@ public:
     virtual Duct* copy() const;
     const Geom& geom() const;
     
+    const pb::Duct& getDuctPb() const { return _ductpb;}
+
     // Solving
     virtual void residual(const TaSystem&,arma::subview_col<d>&& residual) const;
 
diff --git a/src/duct/geom.cpp b/src/duct/geom.cpp
index c580f75..48ff55b 100644
--- a/src/duct/geom.cpp
+++ b/src/duct/geom.cpp
@@ -8,12 +8,14 @@
 #include "geom.h"
 #include "grid.h"
 #include "duct.pb.h"
-#include "tasmet_pyeval.h"
+#include "tasmet_evalscript.h"
 #include <memory>
 #include "tasmet_exception.h"
+#include "tasmet_tracer.h"
 
 Geom::Geom(const pb::Duct& duct) {
 
+    TRACE(15,"Geom::Geom");
     // Step one: create the grid
     std::unique_ptr<Grid> grid;
     switch (duct.gridtype()) {
@@ -35,28 +37,32 @@ Geom::Geom(const pb::Duct& duct) {
     x = grid->getx();
 
     const char* invS = "Invalid cross-sectional area function";
-    EvaluateFun Sfun(duct.area(),invS);
-    Sfun.addGlobalDef("L",duct.length());
-    S = Sfun(x);
-
+    {
+        EvaluateFun Sfun(duct.area(),invS);
+        Sfun.addGlobalDef("L",duct.length());
+        S = Sfun(x);
+    }
     if(min(S) <= 0) {
         throw TaSMETError(invS);
     }
     
 
     const char* invpor = "Invalid porosity function";
-    EvaluateFun phifun(duct.phi(),invpor);
-    Sfun.addGlobalDef("L",duct.length());
-    phi = Sfun(x);
-
+    {
+        EvaluateFun phifun(duct.phi(),invpor);
+        phifun.addGlobalDef("L",duct.length());
+        phi = phifun(x);
+    }
     if(min(phi) <= 0 || max(phi) > 1) {
         throw TaSMETError(invpor);
     }
 
     const char* invrh = "Invalid hydraulic radius function";
-    EvaluateFun rhfun(duct.rh(),invrh);
-    Sfun.addGlobalDef("L",duct.length());
-    rh = Sfun(x);
+    {
+        EvaluateFun rhfun(duct.rh(),invrh);
+        rhfun.addGlobalDef("L",duct.length());
+        rh = rhfun(x);
+    }
     if(min(rh) <= 0 ) {
         throw TaSMETError(invrh);
     }
diff --git a/src/ductbc/adiabaticwall.cpp b/src/ductbc/adiabaticwall.cpp
index 73e2f48..a43a5a2 100644
--- a/src/ductbc/adiabaticwall.cpp
+++ b/src/ductbc/adiabaticwall.cpp
@@ -9,14 +9,14 @@
 #include "ductbc.pb.h"
 #include "tasmet_tracer.h"
 #include "tasmet_assert.h"
-
+#include "tasystem.h"
+#include "duct.h"
 
 AdiabaticWall::AdiabaticWall(const us id,
                        const TaSystem& sys,
                        const pb::DuctBc& dbc):
-    Segment(id,dbc.name()),
-    _side(dbc.side()),
-    _duct_id(dbc.duct_id())
+    DuctBc(id,dbc)
+    
 {
     TRACE(15,"AdiabaticWall(id,sys,dbc)");
     tasmet_assert(dbc.type() == pb::AdiabaticWall,"Wrong type given to constructor");
@@ -24,8 +24,8 @@ AdiabaticWall::AdiabaticWall(const us id,
     
 }
 AdiabaticWall::AdiabaticWall(const AdiabaticWall& o):
-    Segment(o) {
-
+    DuctBc(o)
+{
     TRACE(15,"AdiabaticWall(o)");
 
 }
@@ -39,19 +39,21 @@ vd AdiabaticWall::initialSolution(const TaSystem& sys) const {
     return vd();
 }
 
-void AdiabaticWall::residual(const TaSystem&,
-                          arma::subview_col<d>&& residual
-                          ) const {
+void AdiabaticWall::residual(const TaSystem& sys,
+                             arma::subview_col<d>&& residual
+                             ) const {
 
     TRACE(15,"AdiabaticWall::residual()");
 
-
 }
 
-us AdiabaticWall::getNEqs(const TaSystem&) const {
+us AdiabaticWall::getNEqs(const TaSystem& sys) const {
     TRACE(15,"AdiabaticWall::getNEqs()");    
-
-
+    // u = 0
+    // dT/dx = 0
+    // dTs/dx = 0 => 3 eqs
+    bool has_solideq = getDuct(sys).getDuctPb().stempmodel() != pb::Prescribed;
+    return sys.Ns()*(has_solideq ? 3: 2);
 }
 void AdiabaticWall::show(const TaSystem&,us verbosity_level) const {
     TRACE(15,"AdiabaticWall::show()");    
diff --git a/src/ductbc/adiabaticwall.h b/src/ductbc/adiabaticwall.h
index 110c0f0..63bd04f 100644
--- a/src/ductbc/adiabaticwall.h
+++ b/src/ductbc/adiabaticwall.h
@@ -9,7 +9,7 @@
 #ifndef ADIABATICWALL_H
 #define ADIABATICWALL_H
 #include "segment.h"
-#include "ductbc.pb.h"
+#include "ductbc.h"
 
 class TaSystem;
 class Variable;
@@ -19,7 +19,7 @@ class Variable;
  * blocks the flow and does not allow any axial heat conduction.
  */
 
-class AdiabaticWall: public Segment {
+class AdiabaticWall: public DuctBc {
     pb::DuctSide _side;         /**< Duct side at which this b.c. works */
     us _duct_id;                /**< ID of Duct for this b.c. */
 
diff --git a/src/ductbc/ductbc.cpp b/src/ductbc/ductbc.cpp
index f5266b5..5ebd4fa 100644
--- a/src/ductbc/ductbc.cpp
+++ b/src/ductbc/ductbc.cpp
@@ -11,8 +11,16 @@
 #include "tasmet_assert.h"
 #include "pressurebc.h"
 #include "adiabaticwall.h"
+#include "tasystem.h"
+#include "duct.h"
 
-Segment* DuctBc::newDuctBc(const us id,
+const Duct& DuctBc::getDuct(const TaSystem& sys) const {
+
+    return dynamic_cast<const Duct&>(sys.getSegment(_dbc.duct_id()));
+
+}
+
+DuctBc* DuctBc::newDuctBc(const us id,
                           const TaSystem& sys,
                           const pb::DuctBc& dbc) {
 
diff --git a/src/ductbc/ductbc.h b/src/ductbc/ductbc.h
index e89fab0..d700f87 100644
--- a/src/ductbc/ductbc.h
+++ b/src/ductbc/ductbc.h
@@ -9,17 +9,27 @@
 #ifndef DUCTBC_H
 #define DUCTBC_H
 #include "segment.h"
+#include "ductbc.pb.h"
 
-namespace pb{
-    class DuctBc;
-}
-
-class DuctBc {
+class TaSystem;
+class Duct;
 
+class DuctBc :public Segment {
+    pb::DuctBc _dbc;
 public:
-    static Segment* newDuctBc(const us id,
-                             const TaSystem& sys,
-                             const pb::DuctBc&);
+    DuctBc(const us id,
+           const pb::DuctBc& dbc):
+        Segment(id,dbc.name()),
+        _dbc(dbc)    {}
+
+    DuctBc(const DuctBc& o): Segment(o),_dbc(o._dbc) {}
+
+    static DuctBc* newDuctBc(const us id,
+                              const TaSystem& sys,
+                              const pb::DuctBc&);
+
+    const Duct& getDuct(const TaSystem&) const;
+
 };
 
 
diff --git a/src/ductbc/pressurebc.cpp b/src/ductbc/pressurebc.cpp
index 46536cd..710fe46 100644
--- a/src/ductbc/pressurebc.cpp
+++ b/src/ductbc/pressurebc.cpp
@@ -9,20 +9,19 @@
 #include "pressurebc.h"
 #include "ductbc.pb.h"
 #include "tasmet_tracer.h"
+#include "tasystem.h"
+#include "duct.h"
 
 PressureBc::PressureBc(const us id,
                        const TaSystem& sys,
                        const pb::DuctBc& dbc):
-    Segment(id,dbc.name())
+    DuctBc(id,dbc)
 {
     TRACE(15,"PressureBc(id,sys,dbc)");
 
-
-    
-
 }
 PressureBc::PressureBc(const PressureBc& o):
-    Segment(o) {
+    DuctBc(o) {
 
     TRACE(15,"PressureBc(o)");
 
@@ -46,10 +45,14 @@ void PressureBc::residual(const TaSystem&,
 
 }
 
-us PressureBc::getNEqs(const TaSystem&) const {
+us PressureBc::getNEqs(const TaSystem& sys) const {
     TRACE(15,"PressureBc::getNEqs()");    
-
-
+    // p = x
+    // T = x
+    // This one only if the duct solves for solid
+    // Ts = x => 3 equations
+    bool has_solideq = getDuct(sys).getDuctPb().stempmodel() != pb::Prescribed;
+    return sys.Ns()*(has_solideq ? 3: 2);
 }
 void PressureBc::show(const TaSystem&,us verbosity_level) const {
     TRACE(15,"PressureBc::show()");    
diff --git a/src/ductbc/pressurebc.h b/src/ductbc/pressurebc.h
index 4dccb2f..7a16d7f 100644
--- a/src/ductbc/pressurebc.h
+++ b/src/ductbc/pressurebc.h
@@ -8,14 +8,15 @@
 #pragma once
 #ifndef PRESSUREBC_H
 #define PRESSUREBC_H
-#include "segment.h"
+#include "ductbc.h"
 #include "ductbc.pb.h"
+
 class TaSystem;
 class Variable;
 
-class PressureBc: public Segment {
+class PressureBc: public DuctBc {
     Variable *_p,*_T,*_Ts;
-    us _duct_id;                /**< ID of Duct for this b.c. */
+
     pb::DuctSide _side;         /**< Duct side at which this b.c. works */
 protected:
     PressureBc(const PressureBc&);
diff --git a/src/gui/mainwindow.ui b/src/gui/mainwindow.ui
index d5ef1b1..319def5 100644
--- a/src/gui/mainwindow.ui
+++ b/src/gui/mainwindow.ui
@@ -1,5 +1,6 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <ui version="4.0">
+ <author>J.A. de Jong</author>
  <class>MainWindow</class>
  <widget class="QMainWindow" name="MainWindow">
   <property name="geometry">
@@ -7,7 +8,7 @@
     <x>0</x>
     <y>0</y>
     <width>683</width>
-    <height>488</height>
+    <height>554</height>
    </rect>
   </property>
   <property name="windowTitle">
@@ -63,21 +64,30 @@
                <item row="5" column="0">
                 <widget class="QLabel" name="label_7">
                  <property name="text">
-                  <string>System type:</string>
+                  <string>&amp;System type:</string>
+                 </property>
+                 <property name="buddy">
+                  <cstring>systemtype</cstring>
                  </property>
                 </widget>
                </item>
                <item row="2" column="0">
                 <widget class="QLabel" name="label_5">
                  <property name="text">
-                  <string>Gas type:</string>
+                  <string>&amp;Gas type:</string>
+                 </property>
+                 <property name="buddy">
+                  <cstring>gastype</cstring>
                  </property>
                 </widget>
                </item>
                <item row="3" column="0">
                 <widget class="QLabel" name="label_15">
                  <property name="text">
-                  <string>Reference temperature</string>
+                  <string>&amp;Reference temperature</string>
+                 </property>
+                 <property name="buddy">
+                  <cstring>T0</cstring>
                  </property>
                 </widget>
                </item>
@@ -114,14 +124,20 @@
                <item row="4" column="0">
                 <widget class="QLabel" name="label_17">
                  <property name="text">
-                  <string>Reference pressure</string>
+                  <string>&amp;Reference pressure</string>
+                 </property>
+                 <property name="buddy">
+                  <cstring>p0</cstring>
                  </property>
                 </widget>
                </item>
                <item row="0" column="0">
                 <widget class="QLabel" name="label">
                  <property name="text">
-                  <string>Number of harmonics: </string>
+                  <string>&amp;Number of harmonics: </string>
+                 </property>
+                 <property name="buddy">
+                  <cstring>nf</cstring>
                  </property>
                 </widget>
                </item>
@@ -135,7 +151,10 @@
                <item row="1" column="0">
                 <widget class="QLabel" name="label_4">
                  <property name="text">
-                  <string>Fundamental frequency: </string>
+                  <string>&amp;Fundamental frequency: </string>
+                 </property>
+                 <property name="buddy">
+                  <cstring>freq</cstring>
                  </property>
                 </widget>
                </item>
@@ -192,7 +211,10 @@
                <item row="0" column="0">
                 <widget class="QLabel" name="label_9">
                  <property name="text">
-                  <string>Segment type</string>
+                  <string>S&amp;egment type</string>
+                 </property>
+                 <property name="buddy">
+                  <cstring>segmenttype</cstring>
                  </property>
                 </widget>
                </item>
@@ -202,7 +224,10 @@
                <item row="1" column="0">
                 <widget class="QLabel" name="label_10">
                  <property name="text">
-                  <string>Segment ID</string>
+                  <string>Segment &amp;ID</string>
+                 </property>
+                 <property name="buddy">
+                  <cstring>segmentid</cstring>
                  </property>
                 </widget>
                </item>
@@ -216,7 +241,10 @@
                <item row="2" column="0">
                 <widget class="QLabel" name="label_8">
                  <property name="text">
-                  <string>Segment name</string>
+                  <string>Segment &amp;name</string>
+                 </property>
+                 <property name="buddy">
+                  <cstring>segmentname</cstring>
                  </property>
                 </widget>
                </item>
@@ -237,7 +265,7 @@
                <item row="3" column="0" colspan="2">
                 <widget class="QPushButton" name="addsegment">
                  <property name="text">
-                  <string>Add/edit segment...</string>
+                  <string>&amp;Add/edit segment...</string>
                  </property>
                 </widget>
                </item>
@@ -320,7 +348,7 @@
      <x>0</x>
      <y>0</y>
      <width>683</width>
-     <height>18</height>
+     <height>21</height>
     </rect>
    </property>
    <widget class="QMenu" name="menuFile">
@@ -432,6 +460,21 @@
    </property>
   </action>
  </widget>
+ <tabstops>
+  <tabstop>nf</tabstop>
+  <tabstop>freq</tabstop>
+  <tabstop>gastype</tabstop>
+  <tabstop>T0</tabstop>
+  <tabstop>p0</tabstop>
+  <tabstop>systemtype</tabstop>
+  <tabstop>segmenttype</tabstop>
+  <tabstop>segmentid</tabstop>
+  <tabstop>segmentname</tabstop>
+  <tabstop>addsegment</tabstop>
+  <tabstop>removesegment</tabstop>
+  <tabstop>backlog</tabstop>
+  <tabstop>segoverview</tabstop>
+ </tabstops>
  <resources>
   <include location="tasmet_resources.qrc"/>
  </resources>
diff --git a/src/main.cpp b/src/main.cpp
index 9e84c06..5334189 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -9,10 +9,6 @@
 #include <QCommandLineParser>
 #include <QCommandLineOption>
 
-// For using python from within Qt
-#include <PythonQt.h>
-
-
 #include "tasmet_config.h"
 #include "tasmet_tracer.h"
 #include "gui/mainwindow.h"
@@ -45,18 +41,6 @@ int main(int argc, char *argv[]) {
     // Q_INIT_RESOURCE(application);
     INITTRACE(15);
 
-    // Initialize PythonQt
-    // PythonQt::init(PythonQt::IgnoreSiteModule | PythonQt::RedirectStdOut);
-    PythonQt::init();
-
-    PythonQt* pyqt = PythonQt::self();
-    PythonQtObjectPtr context = pyqt->getMainModule();
-
-    QVariant rv = context.evalScript("from math import *\n");
-    if(pyqt->hadError()) {
-        return -1;
-    }
-    
     QApplication app(argc, argv);
     
     catchUnixSignals({SIGQUIT, SIGINT, SIGTERM, SIGHUP});
@@ -74,8 +58,6 @@ int main(int argc, char *argv[]) {
     // if (!parser.positionalArguments().isEmpty())
     //     mainWin.loadFile(parser.positionalArguments().first());
 
-
-
     TaSMETMainWindow win;
 
     win.show();
diff --git a/src/protobuf/message_tools.cpp b/src/protobuf/message_tools.cpp
index 2265ca9..9cdd97d 100644
--- a/src/protobuf/message_tools.cpp
+++ b/src/protobuf/message_tools.cpp
@@ -12,9 +12,9 @@
 #include "tasmet_exception.h"
 #include "model.pb.h"
 
-#include <google/protobuf/text_format.h>
-
-using google::protobuf::TextFormat;
+#include <google/protobuf/util/message_differencer.h>
+// #include <google/protobuf/text_format.h>
+// using google::protobuf::TextFormat;
 
 template<typename T>
 T loadMessage(const string& filepath) {
@@ -32,11 +32,11 @@ T loadMessage(const string& filepath) {
 
     T sys;
 
-    std::stringstream strStream;
-    strStream << myfile.rdbuf(); //read the file
-    string data = strStream.str();//str holds the content of the file        
+    // std::stringstream strStream;
+    // strStream << myfile.rdbuf(); //read the file
+    // string data = strStream.str();//str holds the content of the file        
         
-    if(!TextFormat::ParseFromString(data,&sys)) {
+    if(!sys.ParseFromIstream(&myfile)) {
         string error = "Invalid TaSMET Model file: ";
         error += filepath;
         throw TaSMETError(error);
@@ -53,27 +53,20 @@ void saveMessage(const string& filepath,const T& sys) {
     if(!sfile.good()){
         throw TaSMETError("Could not write to file");
     }
-    string data;
-    if(TextFormat::PrintToString(sys,&data)) {
-        // Can maybe assign to itself. Which is no problem
-        // according to C++ docs
+
+    // if(TextFormat::PrintToString(sys,&data)) {
+    if(sys.SerializeToOstream(&sfile)) {
         TRACE(15,"Saving file succeeded");
-        sfile << data << std::flush;
-        
-        // Close file here, such that in can be opened to compare
-        // whether the file is still dirty 
-        sfile.close();
     }
     else {
         throw TaSMETError("Could not save model to file");
-        
     }
 }
 
 // // Returns true when the two messages are equal
 template<typename T>
 bool compareMessage(const T& s1,const T& s2) {
-    return (s1.SerializeAsString()==s2.SerializeAsString());
+    return google::protobuf::util::MessageDifferencer::Equals(s1,s2);
 }
 
 // Explicit instantiation for pb::Model
diff --git a/src/protobuf/model.proto b/src/protobuf/model.proto
index 0766be7..b8d24f5 100644
--- a/src/protobuf/model.proto
+++ b/src/protobuf/model.proto
@@ -8,11 +8,7 @@ message Model {
 
     optional System system = 1;
     optional SolverParams sparams = 2;
-
     optional string backlog = 3 [default="Type your info here..."];
-
     repeated double solution = 4;
 
-
-
 }
diff --git a/src/solver/newton_raphson.h b/src/solver/newton_raphson.h
index d55e595..85ac1a0 100644
--- a/src/solver/newton_raphson.h
+++ b/src/solver/newton_raphson.h
@@ -10,16 +10,18 @@
 #define NEWTON_RAPHSON_H
 #include "solver.h"
 
+/**
+ * Newton Raphson Solver implementation.
+ * 
+ */
 class NewtonRaphson: public Solver<GradientNonlinearSystem,vd> {
 
-    d _dampfac = 1.0;
-    d _maxiter = 10000;
+    d _dampfac = 1.0; /**< This is the applied damping factor */
 public:
-    NewtonRaphson(const GradientNonlinearSystem&sys,d dampfac,us maxiter):
+    NewtonRaphson(const GradientNonlinearSystem&sys,d dampfac=1.0):
         Solver(sys),
-        _dampfac(dampfac),
-        _maxiter(maxiter){}
-
+        _dampfac(dampfac)
+        {}
 protected:
     void start_implementation(GradientNonlinearSystem& system,progress_callback* callback);
     
diff --git a/src/solver/solver.cpp b/src/solver/solver.cpp
index e8bb695..b5c46f5 100644
--- a/src/solver/solver.cpp
+++ b/src/solver/solver.cpp
@@ -8,13 +8,32 @@
 #include "tasmet_tracer.h"
 #include "solver.h"
 #include "tasmet_exception.h"
-#include <chrono>
+#include "tasmet_io.h"
+
+SolverAction ostream_progress_callback(SolverProgress pg,std::ostream* out, d funtol,d reltol) {
+
+    // Create reference to ostream
+    std::ostream& myout = (out) ? *out : std::cout;
+
+    myout << "Iteration: " << pg.iteration;
+    myout << " , Residual: " << pg.fun_err << " . Solution error: ";
+    myout << pg.rel_err << endl;
+
+    if(pg.fun_err <= funtol && pg.rel_err <= reltol) {
+        myout << "Solution found within required tolerance" << endl;        
+        return Stop;
+    }
+    
+    return Continue;
+}
 
 
-template<typename system_T,typename result_T>
-static void SolverThread(Solver<system_T,result_T>* solver,
-                         system_T* system,
-                         progress_callback* callback);
+progress_callback simple_progress_callback(d funtol,d reltol,std::ostream* out) {
+    using namespace std::placeholders;  // for _1, _2, _3...
+    using std::bind;
+    return bind(ostream_progress_callback,_1,out,funtol,reltol);
+    
+}
 
 
 template<typename system_T,typename result_T>
@@ -33,28 +52,15 @@ Solver<system_T,result_T>::~Solver(){
 template<typename system_T,typename result_T>
 void Solver<system_T,result_T>::start(progress_callback* callback){
 
-        TRACE(15,"Waiting for solver...");
         start_implementation(*_sys,callback);
 
 }
 
-
 template<typename system_T,typename result_T>
 result_T Solver<system_T,result_T>::getSolution() {
 
     return _sys->getSolution();
 }
-template<typename system_T,typename result_T>
-void SolverThread(Solver<system_T,result_T>* solver,system_T* system,progress_callback* callback) {
-    assert(system);
-    
-    solver->_running = true;
-
-    solver->start_implementation(*system,callback);
-
-    solver->_running = false;
-
-}
 
 
 // Explicit instantiation for some types of systems and results
diff --git a/src/solver/solver.h b/src/solver/solver.h
index 5c7e314..f459ce9 100644
--- a/src/solver/solver.h
+++ b/src/solver/solver.h
@@ -13,44 +13,88 @@
 #include <functional>
 #include <thread>
 #include <atomic>
-
+#include <ostream>
 #include "tasmet_types.h"
 #include "tasmet_assert.h"
 #include "system.h"
 
+/**
+ * Passed to the progress_callback function to keep track of the
+ * Solver progress.
+ * 
+ */
 struct SolverProgress
 {
-    size_t iteration = 0;
+    size_t iteration = 0;       /**< The current iteration no. */
     d fun_err = 1e0;
     d rel_err = 1e0;
-    bool done = false;
+    bool error = false;          /**< Tells the callback that the
+                                   Solver stops as an action of
+                                   itself. Probably due to an internal
+                                   error. */
 };
 
 
+/**
+ * This enum should be return by the Solver callback function.
+ * 
+ */
 enum SolverAction{
-    Continue=0,
-    Stop=1
+    Continue=0,                 /**< Tells the Solver to continue */
+    Stop=1                      /**< Tells the Solver to stop. */
 };
 
+namespace pb{
+    class SolverParams;
+}
+
+/// Typedef for the used callback which is called by a Solver instance.
 typedef std::function<SolverAction(SolverProgress)> progress_callback;
 
+/*! 
+ * Returns a progress calback function that prints the progress
+ * to the given ostream. This is a simple CLI callback function which
+ * prints the solver progress to the given ostream as input. If no
+ * ostream is given, the callback will print to stdout.
+ *
+ * @param funtol The required residual tolerance to stop the solver
+ * @param reltol The required solution tolerance to stop the solver
+ * @param out The ostream to output the progress to.
+ *
+ * @return A std::function object that can be passed to Solver::Start().
+ */
+progress_callback simple_progress_callback(d funtol,d reltol,
+                                           std::ostream* out=nullptr);
+
+
 template<typename system_T,typename result_T>
 class Solver {
 
 protected:
-    system_T* _sys = nullptr;
+    system_T* _sys = nullptr;   /**< Stores and owns a system (creates
+                                   a copy from the given System in the
+                                   constructor */
 public:
     Solver(const system_T& sys);
     Solver(const Solver&)=delete;
     Solver& operator=(const Solver&)=delete;
     
+    /// Start the solver, using the given progress callback
     void start(progress_callback* callback);
 
-    // Returns the solution of the problem
+    /// Returns the solution of the problem after the Solver is done
+    /// solving.
     result_T getSolution();
 
     virtual ~Solver();
+
+    /// Create a new Solver instance based on protobuf
+    /// parameters. Throws a TaSMETError on error.
+    static Solver* newSolver(const pb::SolverParams& params);
+
+
 protected:
+    /// This member fcn should be implemented by the Solver instance.
     virtual void start_implementation(system_T& sys,progress_callback*)=0;
 };
 
diff --git a/src/sys/tasystem.cpp b/src/sys/tasystem.cpp
index a0cfcb8..3f38a44 100644
--- a/src/sys/tasystem.cpp
+++ b/src/sys/tasystem.cpp
@@ -129,8 +129,10 @@ vd TaSystem::residual() const {
     us total_neqs = arma::sum(neqs);
 
     if(total_neqs>constants::maxndofs)      {
-        throw TaSMETError("Too many DOFS required."
-                          " Problem too large.");
+        stringstream error;
+        error << "Too many DOFS required. Problem too large. Number of equations computed: ";
+        error << total_neqs;
+        throw TaSMETError(error);
     }
 
     // This vector of indices stores the last equation number + 1 for
@@ -237,11 +239,12 @@ vus TaSystem::getNDofs() const  {
     return Ndofs;
 }
 vus TaSystem::getNEqs() const  {
-    TRACE(0,"TaSystem::getNDofs()");
+    TRACE(15,"TaSystem::getNEqs()");
     vus Neqs(_segs.size());
     us i=0;
     for (auto seg :_segs) {
         Neqs(i)=seg.second->getNEqs(*this);
+        VARTRACE(15,Neqs(i));
         i++;
     }
     return Neqs;
diff --git a/src/tasmet_evalscript.cpp b/src/tasmet_evalscript.cpp
new file mode 100644
index 0000000..53e8cda
--- /dev/null
+++ b/src/tasmet_evalscript.cpp
@@ -0,0 +1,89 @@
+// tasmet_evalscript.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+
+#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;
+
+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");
+    }
+
+}
+
+EvaluateFun::EvaluateFun(const string& fun_return,
+                         const string& err_msg):
+    _err_msg(err_msg),
+    _fun_return(fun_return)
+{
+    TRACE(15,"EvaluateFun::EvaluateFun()");
+    _chai = getChaiScriptInstance();
+    if(!_chai) throw TaSMETBadAlloc();
+
+    string script = "def myfun(x) {\n";
+    script += "return " + fun_return + "; }\n";
+    
+    wrap_eval<void>(_chai.get(),script);
+    
+}
+void EvaluateFun::addGlobalDef(const string& name,const d value) {
+    TRACE(15,"EvaluateFun::addGlobalDef()");
+
+    _chai->add_global(chaiscript::var(value),name);
+}
+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());
+        
+    }
+    return y;
+}
+
+EvaluateFun::~EvaluateFun() {
+
+}
+
+
+//////////////////////////////////////////////////////////////////////
diff --git a/src/tasmet_evalscript.h b/src/tasmet_evalscript.h
new file mode 100644
index 0000000..77667ed
--- /dev/null
+++ b/src/tasmet_evalscript.h
@@ -0,0 +1,41 @@
+// tasmet_evalscript.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+//
+// Description: This file implements a helper class to evaluate simple
+// 1D single-parameter python math functions
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef TASMET_EVALSCRIPT_H
+#define TASMET_EVALSCRIPT_H
+#include "tasmet_types.h"
+#include <memory>
+
+namespace chaiscript {
+    class ChaiScript;
+}
+
+// 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;
+public:
+    EvaluateFun(const string& fun_return,const string& err_msg = "Script error");
+    
+    // Add a global definition to the namespace
+    void addGlobalDef(const string& name,
+                      const d value);
+    
+    // Evaluate a single function at multiple points and return
+    // the result for each point
+    vd operator()(const vd& points);
+
+    // Used to cleanup the python namespace
+    ~EvaluateFun();
+};
+
+#endif // TASMET_EVALSCRIPT_H
+//////////////////////////////////////////////////////////////////////
diff --git a/src/tasmet_io.h b/src/tasmet_io.h
index 96181ad..cbfc9c0 100644
--- a/src/tasmet_io.h
+++ b/src/tasmet_io.h
@@ -10,6 +10,7 @@
 #define TASMET_IO_H
 #include <iostream>
 #include <iomanip>
+#include <sstream>
 
 //Spoiling global namespace with often used functions and variables 
 using std::cout; /* Output to stdout */				
@@ -20,6 +21,7 @@ using std::cin;		/* Input from stdin */
 using std::ostream; 
 using std::cerr; 
 using std::ios; 
+using std::stringstream;
 
 #endif // TASMET_IO_H
 //////////////////////////////////////////////////////////////////////
diff --git a/src/tasmet_solvemodel.cpp b/src/tasmet_solvemodel.cpp
index c079be5..2d75ddf 100644
--- a/src/tasmet_solvemodel.cpp
+++ b/src/tasmet_solvemodel.cpp
@@ -9,9 +9,10 @@
 #include "tasmet_io.h"
 #include "tasmet_exception.h"
 #include "sys/tasystem.h"
+#include "solver/newton_raphson.h"
 #include "protobuf/model.pb.h"
-// For using python from within Qt
-#include <PythonQt.h>
+#include <memory>
+
 
 void usage(const char* fn) {
     cout << "Usage: " << endl;
@@ -22,18 +23,6 @@ int main(int argc, char *argv[]) {
 
     INITTRACE(10);
 
-    // Initialize PythonQt
-    // PythonQt::init(PythonQt::IgnoreSiteModule | PythonQt::RedirectStdOut);
-    PythonQt::init();
-    
-    PythonQt* pyqt = PythonQt::self();
-    PythonQtObjectPtr context = pyqt->getMainModule();
-
-    QVariant rv = context.evalScript("from math import *\n");
-    if(pyqt->hadError()) {
-        return -1;
-    }
-
     if(argc != 2) {
         usage(argv[0]);
         return -1;
@@ -53,10 +42,10 @@ int main(int argc, char *argv[]) {
     cout << "Loaded model: " << endl;
     cout << model.system().DebugString();
 
-    TaSystem* system;
+    std::unique_ptr<TaSystem> system;
 
     try {
-        system = new TaSystem(model.system());
+        system = std::unique_ptr<TaSystem>(new TaSystem(model.system()));
     }
     catch(TaSMETError &e) {
         cerr << "Model initialization error: " << endl;
@@ -64,6 +53,11 @@ int main(int argc, char *argv[]) {
         return -1;
     }
 
+    NewtonRaphson solver(*system.get());
+
+    progress_callback cb = simple_progress_callback(1e-6,1e-6);
+
+    solver.start(&cb);
 }