diff --git a/apps/sequence/sequence.cpp b/apps/sequence/sequence.cpp
index e9a317f04..4824baf1c 100644
--- a/apps/sequence/sequence.cpp
+++ b/apps/sequence/sequence.cpp
@@ -336,9 +336,7 @@ T Sequence::approximateToNextRank(int n, SequenceContext * sqctx) const {
     {
       ctx.setValueForSymbol(un, &unSymbol);
       ctx.setValueForSymbol(vn, &vnSymbol);
-      Poincare::Complex<T> e = Poincare::Complex<T>::Float(n);
-      ctx.setExpressionForSymbolName(&e, &nSymbol, *sqctx);
-      return expression(sqctx)->template approximateToScalar<T>(ctx);
+      return expression(sqctx)->approximateWithValueForSymbol(symbol(), (T)n, ctx);
     }
     case Type::SingleRecurrence:
     {
@@ -349,9 +347,7 @@ T Sequence::approximateToNextRank(int n, SequenceContext * sqctx) const {
       ctx.setValueForSymbol(unm1, &unSymbol);
       ctx.setValueForSymbol(vn, &vn1Symbol);
       ctx.setValueForSymbol(vnm1, &vnSymbol);
-      Poincare::Complex<T> e = Poincare::Complex<T>::Float(n-1);
-      ctx.setExpressionForSymbolName(&e, &nSymbol, *sqctx);
-      return expression(sqctx)->template approximateToScalar<T>(ctx);
+      return expression(sqctx)->approximateWithValueForSymbol(symbol(), (T)(n-1), ctx);
     }
     default:
     {
@@ -365,9 +361,7 @@ T Sequence::approximateToNextRank(int n, SequenceContext * sqctx) const {
       ctx.setValueForSymbol(unm2, &unSymbol);
       ctx.setValueForSymbol(vnm1, &vn1Symbol);
       ctx.setValueForSymbol(vnm2, &vnSymbol);
-      Poincare::Complex<T> e = Poincare::Complex<T>::Float(n-2);
-      ctx.setExpressionForSymbolName(&e, &nSymbol, *sqctx);
-      return expression(sqctx)->template approximateToScalar<T>(ctx);
+      return expression(sqctx)->approximateWithValueForSymbol(symbol(), (T)(n-2), ctx);
     }
   }
 }
diff --git a/apps/shared/function.cpp b/apps/shared/function.cpp
index 5727f0bea..584f83e03 100644
--- a/apps/shared/function.cpp
+++ b/apps/shared/function.cpp
@@ -41,11 +41,7 @@ void Function::setActive(bool active) {
 
 template<typename T>
 T Function::templatedApproximateAtAbscissa(T x, Poincare::Context * context) const {
-  Poincare::VariableContext<T> variableContext = Poincare::VariableContext<T>(symbol(), context);
-  Poincare::Symbol xSymbol(symbol());
-  Poincare::Complex<T> e = Poincare::Complex<T>::Float(x);
-  variableContext.setExpressionForSymbolName(&e, &xSymbol, variableContext);
-  return expression(context)->approximateToScalar<T>(variableContext);
+  return expression(context)->approximateWithValueForSymbol(symbol(), x, *context);
 }
 
 }
diff --git a/poincare/include/poincare/derivative.h b/poincare/include/poincare/derivative.h
index a24bfd1f4..85e10c3a4 100644
--- a/poincare/include/poincare/derivative.h
+++ b/poincare/include/poincare/derivative.h
@@ -28,8 +28,8 @@ private:
   Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<float>(context, angleUnit); }
   Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const override { return templatedApproximate<double>(context, angleUnit); }
   template<typename T> Expression * templatedApproximate(Context& context, AngleUnit angleUnit) const;
-  template<typename T> T growthRateAroundAbscissa(T x, T h, VariableContext<T> variableContext, AngleUnit angleUnit) const;
-  template<typename T> T riddersApproximation(VariableContext<T> variableContext, AngleUnit angleUnit, T x, T h, T * error) const;
+  template<typename T> T growthRateAroundAbscissa(T x, T h, Context & context, AngleUnit angleUnit) const;
+  template<typename T> T riddersApproximation(Context & context, AngleUnit angleUnit, T x, T h, T * error) const;
   // TODO: Change coefficients?
   constexpr static double k_maxErrorRateOnApproximation = 0.001;
   constexpr static double k_minInitialRate = 0.01;
diff --git a/poincare/include/poincare/expression.h b/poincare/include/poincare/expression.h
index b6252cceb..2bb8e8eeb 100644
--- a/poincare/include/poincare/expression.h
+++ b/poincare/include/poincare/expression.h
@@ -261,6 +261,17 @@ public:
   template<typename T> Expression * approximate(Context& context, AngleUnit angleUnit = AngleUnit::Default) const;
   template<typename T> T approximateToScalar(Context& context, AngleUnit angleUnit = AngleUnit::Default) const;
   template<typename T> static T approximateToScalar(const char * text, Context& context, AngleUnit angleUnit = AngleUnit::Default);
+  template<typename T> T approximateWithValueForSymbol(char symbol, T x, Context & context) const;
+
+  /* Expression roots/extrema solver*/
+  struct Point {
+    double abscissa;
+    double value;
+  };
+  Point nextMinimum(char symbol, double start, double step, double max, Context & context) const;
+  Point nextMaximum(char symbol, double start, double step, double max, Context & context) const;
+  double nextRoot(char symbol, double start, double step, double max, Context & context) const;
+  Point nextIntersection(char symbol, double start, double step, double max, Context & context, const Expression * expression) const;
 protected:
   /* Constructor */
   Expression() : m_parent(nullptr) {}
@@ -318,6 +329,18 @@ private:
   virtual Expression * privateApproximate(SinglePrecision p, Context& context, AngleUnit angleUnit) const = 0;
   virtual Expression * privateApproximate(DoublePrecision p, Context& context, AngleUnit angleUnit) const = 0;
 
+  /* Expression roots/extrema solver*/
+  constexpr static double k_solverPrecision = 1.0E-5;
+  constexpr static double k_sqrtEps = 1.4901161193847656E-8; // sqrt(DBL_EPSILON)
+  constexpr static double k_goldenRatio = 0.381966011250105151795413165634361882279690820194237137864; // (3-sqrt(5))/2
+  typedef double (*EvaluationAtAbscissa)(char symbol, double abscissa, Context & context, const Expression * expression0, const Expression * expression1);
+  Point nextMinimumOfExpression(char symbol, double start, double step, double max, EvaluationAtAbscissa evaluation, Context & context, const Expression * expression = nullptr, bool lookForRootMinimum = false) const;
+  void bracketMinimum(char symbol, double start, double step, double max, double result[3], EvaluationAtAbscissa evaluation, Context & context, const Expression * expression = nullptr) const;
+  Point brentMinimum(char symbol, double ax, double bx, EvaluationAtAbscissa evaluation, Context & context, const Expression * expression = nullptr) const;
+  double nextIntersectionWithExpression(char symbol, double start, double step, double max, EvaluationAtAbscissa evaluation, Context & context, const Expression * expression) const;
+  void bracketRoot(char symbol, double start, double step, double max, double result[2], EvaluationAtAbscissa evaluation, Context & context, const Expression * expression) const;
+  double brentRoot(char symbol, double ax, double bx, double precision, EvaluationAtAbscissa evaluation, Context & context, const Expression * expression) const;
+
   Expression * m_parent;
 };
 
diff --git a/poincare/include/poincare/integral.h b/poincare/include/poincare/integral.h
index c0acca9ce..2c3ab3f98 100644
--- a/poincare/include/poincare/integral.h
+++ b/poincare/include/poincare/integral.h
@@ -34,12 +34,12 @@ private:
   };
   constexpr static int k_maxNumberOfIterations = 10;
 #ifdef LAGRANGE_METHOD
-  template<typename T> T lagrangeGaussQuadrature(T a, T b, VariableContext<T> xContext, AngleUnit angleUnit) const;
+  template<typename T> T lagrangeGaussQuadrature(T a, T b, Context & context, AngleUnit angleUnit) const;
 #else
-  template<typename T> DetailedResult<T> kronrodGaussQuadrature(T a, T b, VariableContext<T> xContext, AngleUnit angleUnit) const;
-  template<typename T> T adaptiveQuadrature(T a, T b, T eps, int numberOfIterations, VariableContext<T> xContext, AngleUnit angleUnit) const;
+  template<typename T> DetailedResult<T> kronrodGaussQuadrature(T a, T b, Context & context, AngleUnit angleUnit) const;
+  template<typename T> T adaptiveQuadrature(T a, T b, T eps, int numberOfIterations, Context & context, AngleUnit angleUnit) const;
 #endif
-  template<typename T> T functionValueAtAbscissa(T x, VariableContext<T> xcontext, AngleUnit angleUnit) const;
+  template<typename T> T functionValueAtAbscissa(T x, Context & xcontext, AngleUnit angleUnit) const;
 };
 
 }
diff --git a/poincare/src/derivative.cpp b/poincare/src/derivative.cpp
index 69103758c..b7cacecdc 100644
--- a/poincare/src/derivative.cpp
+++ b/poincare/src/derivative.cpp
@@ -48,16 +48,9 @@ template<typename T>
 Expression * Derivative::templatedApproximate(Context& context, AngleUnit angleUnit) const {
   static T min = sizeof(T) == sizeof(double) ? DBL_MIN : FLT_MIN;
   static T epsilon = sizeof(T) == sizeof(double) ? DBL_EPSILON : FLT_EPSILON;
-  VariableContext<T> xContext = VariableContext<T>('x', &context);
   Symbol xSymbol('x');
-  Expression * xInput = operand(1)->approximate<T>(context, angleUnit);
-  T x = xInput->type() == Type::Complex ? static_cast<Complex<T> *>(xInput)->toScalar() : NAN;
-  delete xInput;
-  Complex<T> e = Complex<T>::Float(x);
-  xContext.setExpressionForSymbolName(&e, &xSymbol, xContext);
-  Expression * fInput = operand(0)->approximate<T>(xContext, angleUnit);
-  T functionValue = fInput->type() == Type::Complex ? static_cast<Complex<T> *>(fInput)->toScalar() : NAN;
-  delete fInput;
+  T x = operand(1)->approximateToScalar<T>(context, angleUnit);
+  T functionValue = operand(0)->approximateWithValueForSymbol('x', x, context);
 
   // No complex/matrix version of Derivative
   if (std::isnan(x) || std::isnan(functionValue)) {
@@ -67,7 +60,7 @@ Expression * Derivative::templatedApproximate(Context& context, AngleUnit angleU
   T error, result;
   T h = k_minInitialRate;
   do {
-    result = riddersApproximation(xContext, angleUnit, x, h, &error);
+    result = riddersApproximation(context, angleUnit, x, h, &error);
     h /= 10.0;
   } while ((std::fabs(error/result) > k_maxErrorRateOnApproximation || std::isnan(error)) && h >= epsilon);
 
@@ -83,23 +76,14 @@ Expression * Derivative::templatedApproximate(Context& context, AngleUnit angleU
 }
 
 template<typename T>
-T Derivative::growthRateAroundAbscissa(T x, T h, VariableContext<T> xContext, AngleUnit angleUnit) const {
-  Symbol xSymbol('x');
-  Complex<T> e = Complex<T>::Float(x + h);
-  xContext.setExpressionForSymbolName(&e, &xSymbol, xContext);
-  Expression * fInput = operand(0)->approximate<T>(xContext, angleUnit);
-  T expressionPlus = fInput->type() == Type::Complex ? static_cast<Complex<T> *>(fInput)->toScalar() : NAN;
-  delete fInput;
-  e = Complex<T>::Float(x-h);
-  xContext.setExpressionForSymbolName(&e, &xSymbol, xContext);
-  fInput = operand(0)->approximate<T>(xContext, angleUnit);
-  T expressionMinus = fInput->type() == Type::Complex ? static_cast<Complex<T> *>(fInput)->toScalar() : NAN;
-  delete fInput;
+T Derivative::growthRateAroundAbscissa(T x, T h, Context & context, AngleUnit angleUnit) const {
+  T expressionPlus = operand(0)->approximateWithValueForSymbol('x', x+h, context);
+  T expressionMinus = operand(0)->approximateWithValueForSymbol('x', x-h, context);
   return (expressionPlus - expressionMinus)/(2*h);
 }
 
 template<typename T>
-T Derivative::riddersApproximation(VariableContext<T> xContext, AngleUnit angleUnit, T x, T h, T * error) const {
+T Derivative::riddersApproximation(Context & context, AngleUnit angleUnit, T x, T h, T * error) const {
   /* Ridders' Algorithm
    * Blibliography:
    * - Ridders, C.J.F. 1982, Advances in Engineering Software, vol. 4, no. 2,
@@ -119,7 +103,7 @@ T Derivative::riddersApproximation(VariableContext<T> xContext, AngleUnit angleU
       a[i][j] = 1;
     }
   }
-  a[0][0] = growthRateAroundAbscissa(x, hh, xContext, angleUnit);
+  a[0][0] = growthRateAroundAbscissa(x, hh, context, angleUnit);
   T ans = 0;
   T errt = 0;
   /* Loop on i: change the step size */
@@ -128,7 +112,7 @@ T Derivative::riddersApproximation(VariableContext<T> xContext, AngleUnit angleU
     /* Make hh an exactly representable number */
     volatile T temp =  x+hh;
     hh = temp - x;
-    a[0][i] = growthRateAroundAbscissa(x, hh, xContext, angleUnit);
+    a[0][i] = growthRateAroundAbscissa(x, hh, context, angleUnit);
     T fac = k_rateStepSize*k_rateStepSize;
     /* Loop on j: compute extrapolation for several orders */
     for (int j = 1; j < 10; j++) {
diff --git a/poincare/src/expression.cpp b/poincare/src/expression.cpp
index 388d258fe..9f896c5ae 100644
--- a/poincare/src/expression.cpp
+++ b/poincare/src/expression.cpp
@@ -11,7 +11,9 @@
 #include <poincare/matrix.h>
 #include <poincare/complex.h>
 #include <poincare/opposite.h>
+#include <poincare/variable_context.h>
 #include <cmath>
+#include <float.h>
 #include "expression_parser.hpp"
 #include "expression_lexer.hpp"
 
@@ -456,6 +458,324 @@ template<typename T> T Expression::epsilon() {
   return epsilon;
 }
 
+/* Expression roots/extrema solver*/
+
+Expression::Point Expression::nextMinimum(char symbol, double start, double step, double max, Context & context) const {
+  return nextMinimumOfExpression(symbol, start, step, max, [](char symbol, double x, Context & context, const Expression * expression0, const Expression * expression1 = nullptr) {
+        return expression0->approximateWithValueForSymbol(symbol, x, context);
+      }, context);
+}
+
+Expression::Point Expression::nextMaximum(char symbol, double start, double step, double max, Context & context) const {
+  Point minimumOfOpposite = nextMinimumOfExpression(symbol, start, step, max, [](char symbol, double x, Context & context, const Expression * expression0, const Expression * expression1 = nullptr) {
+        return -expression0->approximateWithValueForSymbol(symbol, x, context);
+      }, context);
+  return {.abscissa = minimumOfOpposite.abscissa, .value = -minimumOfOpposite.value};
+}
+
+double Expression::nextRoot(char symbol, double start, double step, double max, Context & context) const {
+  return nextIntersectionWithExpression(symbol, start, step, max, [](char symbol, double x, Context & context, const Expression * expression0, const Expression * expression1 = nullptr) {
+        return expression0->approximateWithValueForSymbol(symbol, x, context);
+      }, context, nullptr);
+}
+
+Expression::Point Expression::nextIntersection(char symbol, double start, double step, double max, Poincare::Context & context, const Expression * expression) const {
+  double resultAbscissa = nextIntersectionWithExpression(symbol, start, step, max, [](char symbol, double x, Context & context, const Expression * expression0, const Expression * expression1) {
+        return expression0->approximateWithValueForSymbol(symbol, x, context)-expression1->approximateWithValueForSymbol(symbol, x, context);
+      }, context, expression);
+  Expression::Point result = {.abscissa = resultAbscissa, .value = approximateWithValueForSymbol(symbol, resultAbscissa, context)};
+  if (std::fabs(result.value) < step*k_solverPrecision) {
+    result.value = 0.0;
+  }
+  return result;
+}
+
+Expression::Point Expression::nextMinimumOfExpression(char symbol, double start, double step, double max, EvaluationAtAbscissa evaluate, Context & context, const Expression * expression, bool lookForRootMinimum) const {
+  double bracket[3];
+  Point result = {.abscissa = NAN, .value = NAN};
+  double x = start;
+  bool endCondition = false;
+  do {
+    bracketMinimum(symbol, x, step, max, bracket, evaluate, context, expression);
+    result = brentMinimum(symbol, bracket[0], bracket[2], evaluate, context, expression);
+    x = bracket[1];
+    endCondition = std::isnan(result.abscissa) && (step > 0.0 ? x <= max : x >= max);
+    if (lookForRootMinimum) {
+      endCondition |= std::fabs(result.value) >= k_sqrtEps && (step > 0.0 ? x <= max : x >= max);
+    }
+  } while (endCondition);
+
+  if (std::fabs(result.abscissa) < step*k_solverPrecision) {
+    result.abscissa = 0;
+    result.value = evaluate(symbol, 0, context, this, expression);
+  }
+  if (std::fabs(result.value) < step*k_solverPrecision) {
+    result.value = 0;
+  }
+  if (lookForRootMinimum) {
+    result.abscissa = std::fabs(result.value) >= k_sqrtEps ? NAN : result.abscissa;
+  }
+  return result;
+}
+
+void Expression::bracketMinimum(char symbol, double start, double step, double max, double result[3], EvaluationAtAbscissa evaluate, Context & context, const Expression * expression) const {
+  Point p[3];
+  p[0] = {.abscissa = start, .value = evaluate(symbol, start, context, this, expression)};
+  p[1] = {.abscissa = start+step, .value = evaluate(symbol, start+step, context, this, expression)};
+  double x = start+2.0*step;
+  while (step > 0.0 ? x <= max : x >= max) {
+    p[2] = {.abscissa = x, .value = evaluate(symbol, x, context, this, expression)};
+    if (p[0].value > p[1].value && p[2].value > p[1].value) {
+      result[0] = p[0].abscissa;
+      result[1] = p[1].abscissa;
+      result[2] = p[2].abscissa;
+      return;
+    }
+    if (p[0].value > p[1].value && p[1].value == p[2].value) {
+    } else {
+      p[0] = p[1];
+      p[1] = p[2];
+    }
+    x += step;
+  }
+  result[0] = NAN;
+  result[1] = NAN;
+  result[2] = NAN;
+}
+
+Expression::Point Expression::brentMinimum(char symbol, double ax, double bx, EvaluationAtAbscissa evaluate, Context & context, const Expression * expression) const {
+  /* Bibliography: R. P. Brent, Algorithms for finding zeros and extrema of
+   * functions without calculating derivatives */
+  if (ax > bx) {
+    return brentMinimum(symbol, bx, ax, evaluate, context, expression);
+  }
+  double e = 0.0;
+  double a = ax;
+  double b = bx;
+  double x = a+k_goldenRatio*(b-a);
+  double v = x;
+  double w = x;
+  double fx = evaluate(symbol, x, context, this, expression);
+  double fw = fx;
+  double fv = fw;
+
+  double d = NAN;
+  double u, fu;
+
+  for (int i = 0; i < 100; i++) {
+    double m = 0.5*(a+b);
+    double tol1 = k_sqrtEps*std::fabs(x)+1E-10;
+    double tol2 = 2.0*tol1;
+    if (std::fabs(x-m) <= tol2-0.5*(b-a))  {
+      double middleFax = evaluate(symbol, (x+a)/2.0, context, this, expression);
+      double middleFbx = evaluate(symbol, (x+b)/2.0, context, this, expression);
+      double fa = evaluate(symbol, a, context, this, expression);
+      double fb = evaluate(symbol, b, context, this, expression);
+      if (middleFax-fa <= k_sqrtEps && fx-middleFax <= k_sqrtEps && fx-middleFbx <= k_sqrtEps && middleFbx-fb <= k_sqrtEps) {
+        Point result = {.abscissa = x, .value = fx};
+        return result;
+      }
+    }
+    double p = 0;
+    double q = 0;
+    double r = 0;
+    if (std::fabs(e) > tol1) {
+      r = (x-w)*(fx-fv);
+      q = (x-v)*(fx-fw);
+      p = (x-v)*q -(x-w)*r;
+      q = 2.0*(q-r);
+      if (q>0.0) {
+        p = -p;
+      } else {
+        q = -q;
+      }
+      r = e;
+      e = d;
+    }
+    if (std::fabs(p) < std::fabs(0.5*q*r) && p<q*(a-x) && p<q*(b-x)) {
+      d = p/q;
+      u= x+d;
+      if (u-a < tol2 || b-u < tol2) {
+        d = x < m ? tol1 : -tol1;
+      }
+    } else {
+      e = x<m ? b-x : a-x;
+      d = k_goldenRatio*e;
+    }
+    u = x + (std::fabs(d) >= tol1 ? d : (d>0 ? tol1 : -tol1));
+    fu = evaluate(symbol, u, context, this, expression);
+    if (fu <= fx) {
+      if (u<x) {
+        b = x;
+      } else {
+        a = x;
+      }
+      v = w;
+      fv = fw;
+      w = x;
+      fw = fx;
+      x = u;
+      fx = fu;
+    } else {
+      if (u<x) {
+        a = u;
+      } else {
+        b = u;
+      }
+      if (fu <= fw || w == x) {
+        v = w;
+        fv = fw;
+        w = u;
+        fw = fu;
+      } else if (fu <= fv || v == x || v == w) {
+        v = u;
+        fv = fu;
+      }
+    }
+  }
+  Point result = {.abscissa = NAN, .value = NAN};
+  return result;
+}
+
+double Expression::nextIntersectionWithExpression(char symbol, double start, double step, double max, EvaluationAtAbscissa evaluation, Context & context, const Expression * expression) const {
+  double bracket[2];
+  double result = NAN;
+  static double precisionByGradUnit = 1E6;
+  double x = start+step;
+  do {
+    bracketRoot(symbol, x, step, max, bracket, evaluation, context, expression);
+    result = brentRoot(symbol, bracket[0], bracket[1], std::fabs(step/precisionByGradUnit), evaluation, context, expression);
+    x = bracket[1];
+  } while (std::isnan(result) && (step > 0.0 ? x <= max : x >= max));
+
+  double extremumMax = std::isnan(result) ? max : result;
+  Point resultExtremum[2] = {
+    nextMinimumOfExpression(symbol, start, step, extremumMax, [](char symbol, double x, Context & context, const Expression * expression0, const Expression * expression1) {
+        if (expression1) {
+          return expression0->approximateWithValueForSymbol(symbol, x, context)-expression1->approximateWithValueForSymbol(symbol, x, context);
+        } else {
+          return expression0->approximateWithValueForSymbol(symbol, x, context);
+        }
+      }, context, expression, true),
+    nextMinimumOfExpression(symbol, start, step, extremumMax, [](char symbol, double x, Context & context, const Expression * expression0, const Expression * expression1) {
+        if (expression1) {
+          return expression1->approximateWithValueForSymbol(symbol, x, context)-expression0->approximateWithValueForSymbol(symbol, x, context);
+        } else {
+          return -expression0->approximateWithValueForSymbol(symbol, x, context);
+        }
+      }, context, expression, true)};
+  for (int i = 0; i < 2; i++) {
+    if (!std::isnan(resultExtremum[i].abscissa) && (std::isnan(result) || std::fabs(result - start) > std::fabs(resultExtremum[i].abscissa - start))) {
+      result = resultExtremum[i].abscissa;
+    }
+  }
+  if (std::fabs(result) < step*k_solverPrecision) {
+    result = 0;
+  }
+  return result;
+}
+
+void Expression::bracketRoot(char symbol, double start, double step, double max, double result[2], EvaluationAtAbscissa evaluation, Context & context, const Expression * expression) const {
+  double a = start;
+  double b = start+step;
+  while (step > 0.0 ? b <= max : b >= max) {
+    double fa = evaluation(symbol, a, context, this, expression);
+    double fb = evaluation(symbol, b, context, this, expression);
+    if (fa*fb <= 0) {
+      result[0] = a;
+      result[1] = b;
+      return;
+    }
+    a = b;
+    b = b+step;
+  }
+  result[0] = NAN;
+  result[1] = NAN;
+}
+
+
+double Expression::brentRoot(char symbol, double ax, double bx, double precision, EvaluationAtAbscissa evaluation, Context & context, const Expression * expression) const {
+  if (ax > bx) {
+    return brentRoot(symbol, bx, ax, precision, evaluation, context, expression);
+  }
+  double a = ax;
+  double b = bx;
+  double c = bx;
+  double d = b-a;
+  double e = b-a;
+  double fa = evaluation(symbol, a, context, this, expression);
+  double fb = evaluation(symbol, b, context, this, expression);
+  double fc = fb;
+  for (int i = 0; i < 100; i++) {
+    if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
+      c = a;
+      fc = fa;
+      e = b-a;
+      d = b-a;
+    }
+    if (std::fabs(fc) < std::fabs(fb)) {
+      a = b;
+      b = c;
+      c = a;
+      fa = fb;
+      fb = fc;
+      fc = fa;
+    }
+    double tol1 = 2.0*DBL_EPSILON*std::fabs(b)+0.5*precision;
+    double xm = 0.5*(c-b);
+    if (std::fabs(xm) <= tol1 || fb == 0.0) {
+      double fbcMiddle = evaluation(symbol, 0.5*(b+c), context, this, expression);
+      double isContinuous = (fb <= fbcMiddle && fbcMiddle <= fc) || (fc <= fbcMiddle && fbcMiddle <= fb);
+      if (isContinuous) {
+        return b;
+      }
+    }
+    if (std::fabs(e) >= tol1 && std::fabs(fa) > std::fabs(b)) {
+      double s = fb/fa;
+      double p = 2.0*xm*s;
+      double q = 1.0-s;
+      if (a != c) {
+        q = fa/fc;
+        double r = fb/fc;
+        p = s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
+        q = (q-1.0)*(r-1.0)*(s-1.0);
+      }
+      q = p > 0.0 ? -q : q;
+      p = std::fabs(p);
+      double min1 = 3.0*xm*q-std::fabs(tol1*q);
+      double min2 = std::fabs(e*q);
+      if (2.0*p < (min1 < min2 ? min1 : min2)) {
+        e = d;
+        d = p/q;
+      } else {
+        d = xm;
+        e =d;
+      }
+    } else {
+      d = xm;
+      e = d;
+    }
+    a = b;
+    fa = fb;
+    if (std::fabs(d) > tol1) {
+      b += d;
+    } else {
+      b += xm > 0.0 ? tol1 : tol1;
+    }
+    fb = evaluation(symbol, b, context, this, expression);
+  }
+  return NAN;
+}
+
+template<typename T>
+T Expression::approximateWithValueForSymbol(char symbol, T x, Context & context) const {
+  VariableContext<T> variableContext = VariableContext<T>(symbol, &context);
+  Symbol xSymbol(symbol);
+  Complex<T> e = Complex<T>::Float(x);
+  variableContext.setExpressionForSymbolName(&e, &xSymbol, variableContext);
+  return approximateToScalar<T>(variableContext);
+}
+
 }
 
 template Poincare::Expression * Poincare::Expression::approximate<double>(Context& context, AngleUnit angleUnit) const;
@@ -466,3 +786,5 @@ template double Poincare::Expression::approximateToScalar<double>(Poincare::Cont
 template float Poincare::Expression::approximateToScalar<float>(Poincare::Context&, Poincare::Expression::AngleUnit) const;
 template double Poincare::Expression::epsilon<double>();
 template float Poincare::Expression::epsilon<float>();
+template double Poincare::Expression::approximateWithValueForSymbol(char, double, Poincare::Context&) const;
+template float Poincare::Expression::approximateWithValueForSymbol(char, float, Poincare::Context&) const;
diff --git a/poincare/src/integral.cpp b/poincare/src/integral.cpp
index 8536c9cb5..aa91f3495 100644
--- a/poincare/src/integral.cpp
+++ b/poincare/src/integral.cpp
@@ -49,7 +49,6 @@ Expression * Integral::shallowReduce(Context& context, AngleUnit angleUnit) {
 
 template<typename T>
 Complex<T> * Integral::templatedApproximate(Context & context, AngleUnit angleUnit) const {
-  VariableContext<T> xContext = VariableContext<T>('x', &context);
   Expression * aInput = operand(1)->approximate<T>(context, angleUnit);
   T a = aInput->type() == Type::Complex ? static_cast<Complex<T> *>(aInput)->toScalar() : NAN;
   delete aInput;
@@ -60,9 +59,9 @@ Complex<T> * Integral::templatedApproximate(Context & context, AngleUnit angleUn
     return new Complex<T>(Complex<T>::Float(NAN));
   }
 #ifdef LAGRANGE_METHOD
-  T result = lagrangeGaussQuadrature<T>(a, b, xContext, angleUnit);
+  T result = lagrangeGaussQuadrature<T>(a, b, context, angleUnit);
 #else
-  T result = adaptiveQuadrature<T>(a, b, 0.1, k_maxNumberOfIterations, xContext, angleUnit);
+  T result = adaptiveQuadrature<T>(a, b, 0.1, k_maxNumberOfIterations, context, angleUnit);
 #endif
   return new Complex<T>(Complex<T>::Float(result));
 }
@@ -78,20 +77,14 @@ ExpressionLayout * Integral::privateCreateLayout(PrintFloat::Mode floatDisplayMo
 }
 
 template<typename T>
-T Integral::functionValueAtAbscissa(T x, VariableContext<T> xContext, AngleUnit angleUnit) const {
-  Complex<T> e = Complex<T>::Float(x);
-  Symbol xSymbol('x');
-  xContext.setExpressionForSymbolName(&e, &xSymbol, xContext);
-  Expression * f = operand(0)->approximate<T>(xContext, angleUnit);
-  T result = f->type() == Type::Complex ? static_cast<Complex<T> *>(f)->toScalar() : NAN;
-  delete f;
-  return result;
+T Integral::functionValueAtAbscissa(T x, Context & context, AngleUnit angleUnit) const {
+  return operand(0)->approximateWithValueForSymbol('x', x, context);
 }
 
 #ifdef LAGRANGE_METHOD
 
 template<typename T>
-T Integral::lagrangeGaussQuadrature(T a, T b, VariableContext<T> xContext, AngleUnit angleUnit) const {
+T Integral::lagrangeGaussQuadrature(T a, T b, Context & context, AngleUnit angleUnit) const {
   /* We here use Gauss-Legendre quadrature with n = 5
    * Gauss-Legendre abscissae and weights can be found in
    * C/C++ library source code. */
@@ -105,7 +98,7 @@ T Integral::lagrangeGaussQuadrature(T a, T b, VariableContext<T> xContext, Angle
   T result = 0;
   for (int j = 0; j < 10; j++) {
     T dx = xr * x[j];
-    result += w[j]*(functionValueAtAbscissa(xm+dx, xContext, angleUnit) + functionValueAtAbscissa(xm-dx, xContext, angleUnit));
+    result += w[j]*(functionValueAtAbscissa(xm+dx, context, angleUnit) + functionValueAtAbscissa(xm-dx, context, angleUnit));
   }
   result *= xr;
   return result;
@@ -114,7 +107,7 @@ T Integral::lagrangeGaussQuadrature(T a, T b, VariableContext<T> xContext, Angle
 #else
 
 template<typename T>
-Integral::DetailedResult<T> Integral::kronrodGaussQuadrature(T a, T b, VariableContext<T> xContext, AngleUnit angleUnit) const {
+Integral::DetailedResult<T> Integral::kronrodGaussQuadrature(T a, T b, Context & context, AngleUnit angleUnit) const {
   static T epsilon = sizeof(T) == sizeof(double) ? DBL_EPSILON : FLT_EPSILON;
   static T max = sizeof(T) == sizeof(double) ? DBL_MAX : FLT_MAX;
   /* We here use Kronrod-Legendre quadrature with n = 21
@@ -137,14 +130,14 @@ Integral::DetailedResult<T> Integral::kronrodGaussQuadrature(T a, T b, VariableC
   T dhlgth = std::fabs(hlgth);
 
   T resg = 0;
-  T fc = functionValueAtAbscissa(centr, xContext, angleUnit);
+  T fc = functionValueAtAbscissa(centr, context, angleUnit);
   T resk = wgk[10]*fc;
   T resabs = std::fabs(resk);
   for (int j = 0; j < 5; j++) {
     int jtw = 2*j+1;
     T absc = hlgth*xgk[jtw];
-    T fval1 = functionValueAtAbscissa(centr-absc, xContext, angleUnit);
-    T fval2 = functionValueAtAbscissa(centr+absc, xContext, angleUnit);
+    T fval1 = functionValueAtAbscissa(centr-absc, context, angleUnit);
+    T fval2 = functionValueAtAbscissa(centr+absc, context, angleUnit);
     fv1[jtw] = fval1;
     fv2[jtw] = fval2;
     T fsum = fval1+fval2;
@@ -155,8 +148,8 @@ Integral::DetailedResult<T> Integral::kronrodGaussQuadrature(T a, T b, VariableC
   for (int j = 0; j < 5; j++) {
     int jtwm1 = 2*j;
     T absc = hlgth*xgk[jtwm1];
-    T fval1 = functionValueAtAbscissa(centr-absc, xContext, angleUnit);
-    T fval2 = functionValueAtAbscissa(centr+absc, xContext, angleUnit);
+    T fval1 = functionValueAtAbscissa(centr-absc, context, angleUnit);
+    T fval2 = functionValueAtAbscissa(centr+absc, context, angleUnit);
     fv1[jtwm1] = fval1;
     fv2[jtwm1] = fval2;
     T fsum = fval1+fval2;
@@ -185,17 +178,17 @@ Integral::DetailedResult<T> Integral::kronrodGaussQuadrature(T a, T b, VariableC
 }
 
 template<typename T>
-T Integral::adaptiveQuadrature(T a, T b, T eps, int numberOfIterations, VariableContext<T> xContext, AngleUnit angleUnit) const {
+T Integral::adaptiveQuadrature(T a, T b, T eps, int numberOfIterations, Context & context, AngleUnit angleUnit) const {
   if (shouldStopProcessing()) {
     return NAN;
   }
-  DetailedResult<T> quadKG = kronrodGaussQuadrature(a, b, xContext, angleUnit);
+  DetailedResult<T> quadKG = kronrodGaussQuadrature(a, b, context, angleUnit);
   T result = quadKG.integral;
   if (quadKG.absoluteError <= eps) {
     return result;
   } else if (--numberOfIterations > 0) {
     T m = (a+b)/2;
-    return adaptiveQuadrature<T>(a, m, eps/2, numberOfIterations, xContext, angleUnit) + adaptiveQuadrature<T>(m, b, eps/2, numberOfIterations, xContext, angleUnit);
+    return adaptiveQuadrature<T>(a, m, eps/2, numberOfIterations, context, angleUnit) + adaptiveQuadrature<T>(m, b, eps/2, numberOfIterations, context, angleUnit);
   } else {
     return NAN;
   }