diff --git a/contrib/MathEval/node.cpp b/contrib/MathEval/node.cpp
index 7ced8e2c78be1768e9febc5c8a0366d77984afa3..5d7bc033af63bda830f977a520c46309988b8d5d 100644
--- a/contrib/MathEval/node.cpp
+++ b/contrib/MathEval/node.cpp
@@ -502,6 +502,43 @@ node_derivative(Node * node, char *name, SymbolTable * symbol_table)
      */
     else if (!strcmp(node->data.function.record->name, "Fabs"))
       return node_create('b', '/', node_create('b', '*', node_derivative(node->data.function.child, name, symbol_table), node_copy(node->data.function.child)), node_create('f', symbol_table_lookup(symbol_table, "Sqrt"), node_create('b', '^', node_copy(node->data.function.child), node_create('c', 2.0))));
+
+     /*
+     * Apply rule of step function derivative. 
+     */
+    else if (!strcmp(node->data.function.record->name, "Step"))
+      return node_create('b', '*',node_derivative(node->data.function.child, name, symbol_table), node_create('f', symbol_table_lookup(symbol_table,"Delta"), node_copy(node->data.function.child)));
+
+    /*
+     * Apply rule of floor function derivative. 
+     */
+    else if (!strcmp(node->data.function.record->name, "Floor"))
+       return node_create('b', '*', node_derivative(node->data.function.child, name, symbol_table), node_create('f',symbol_table_lookup(symbol_table, "Delta"), node_copy(node->data.function.child)));
+
+    /*
+     * Apply rule of smoothed heaviside function derivative. 
+     */
+    else if (!strcmp(node->data.function.record->name, "HeavS")){
+      //printf("********************************  Derivative of smoothed heaviside function not implemented yet, return a Delta function ...");
+      return node_create('b', '*',node_derivative(node->data.function.child, name, symbol_table), node_create('f', symbol_table_lookup(symbol_table,"Delta"), node_copy(node->data.function.child)));
+    }
+
+    /*
+     * Apply rule of delta function derivative. 
+     */
+    else if (!strcmp(node->data.function.record->name, "Delta"))
+      return node_create('b', '*',node_derivative(node->data.function.child, name,symbol_table),node_create('f',symbol_table_lookup(symbol_table, "NanDelta"), node_copy(node->data.function.child)));
+     
+
+    /*
+     * Apply rule of nandelta function derivative. 
+     */
+    else if (!strcmp(node->data.function.record->name, "NanDelta"))
+      return node_create('b', '*', node_derivative(node->data.function.child, name, symbol_table), node_create('f',symbol_table_lookup(symbol_table, "NanDelta"), node_copy(node->data.function.child)));
+
+
+
+
     
   case 'u':
     switch (node->data.un_op.operatorr) {
diff --git a/contrib/MathEval/scanner.l b/contrib/MathEval/scanner.l
index 759dad1fe2cf93e900dd2d20979fee920e913e47..69a60418f5d873b43fb1887a2aaf2e2f9ff4819f 100644
--- a/contrib/MathEval/scanner.l
+++ b/contrib/MathEval/scanner.l
@@ -46,7 +46,7 @@ static int input_from_string (char *buffer, int max_size);
 whitespace [\ \t]+
 digit [0-9]
 number ({digit}+|{digit}+"."{digit}*|{digit}*"."{digit}+)([Ee][-+]?{digit}+)?
-function "Exp"|"Log"|"Sqrt"|"Sin"|"Cos"|"Tan"|"Ctan"|"Asin"|"Acos"|"Atan"|"Actan"|"Sinh"|"Cosh"|"Tanh"|"Ctanh"|"Asinh"|"Acosh"|"Atanh"|"Actanh"|"Fabs"|"Rand"|"Log10"
+function "Exp"|"Log"|"Sqrt"|"Sin"|"Cos"|"Tan"|"Ctan"|"Asin"|"Acos"|"Atan"|"Actan"|"Sinh"|"Cosh"|"Tanh"|"Ctanh"|"Asinh"|"Acosh"|"Atanh"|"Actanh"|"Fabs"|"Rand"|"Log10"|"Step"|"Delta"|"NanDelta"|"Floor"|"HeavS"
 identifier [a-zA-Z\_][a-zA-Z0-9\_]*
 
 %%
diff --git a/contrib/MathEval/symbol_table.cpp b/contrib/MathEval/symbol_table.cpp
index dc7db03fcae3fd77849c2ff9edb87e0040a29538..5819936fccb61420c5d533017ac73377b685ed4f 100644
--- a/contrib/MathEval/symbol_table.cpp
+++ b/contrib/MathEval/symbol_table.cpp
@@ -26,6 +26,7 @@
 #include "common.h"
 #include "symbol_table.h"
 #include "xmath.h"
+#include <iostream>
 
 /*
  * Type definition for function accepting single argument of double type and
@@ -43,11 +44,11 @@ symbol_table_create(int length)
   static char    *names[] = {"Exp", "Log", "Sqrt", "Sin", "Cos", "Tan", "Ctan", 
 			     "Asin", "Acos", "Atan", "Actan", "Sinh", "Cosh", "Tanh",
 			     "Ctanh", "Asinh", "Acosh", "Atanh", "Actanh", "Fabs", 
-			     "Rand", "Log10" };
+			     "Rand", "Log10" , "Step", "Delta", "NanDelta", "Floor", "HeavS"};
   static double   (*functions[]) (double) = { exp, log, sqrt, sin, cos, tan, x_ctan, 
 					      asin, acos, atan, x_actan, sinh, cosh, tanh, 
 					      x_ctanh, x_asinh, x_acosh, x_atanh, x_actanh,
-					      fabs, x_rand, log10};
+					      fabs, x_rand, log10, x_step, x_delta , x_nandelta, floor, x_heavs};
   unsigned int i;
   
   /*
diff --git a/contrib/MathEval/xmath.cpp b/contrib/MathEval/xmath.cpp
index 9f431af284e6d372bf164f9368bf1af4817ae41e..5e2cb808ce29afa0a813a67820ab146f8f8cf5a9 100644
--- a/contrib/MathEval/xmath.cpp
+++ b/contrib/MathEval/xmath.cpp
@@ -23,6 +23,10 @@
 #include "xmath.h"
 #include <math.h>
 #include <stdlib.h>
+#include <cmath>
+#include <stdio.h>
+
+using std::fabs;
 
 double
 x_ctan(double x)
@@ -99,3 +103,40 @@ x_rand(double x)
   
   return x*(double)rand()/(double)RAND_MAX;
 }
+double
+x_step(double x)
+{
+	/* 
+	 * Calculate step function value.
+	 */
+	return (x < 0) ? 0 : 1;
+}
+
+double
+x_heavs(double x) //, double eps)
+{
+	/* 
+	 * Calculate smoothed heaviside function value.
+	 */
+  double eps=3.0;
+         return (fabs(x) < eps) ? 0.5*(1+ x/eps +1/3.14*sin(3.14*x/eps)):(x<0? 0:1);
+}
+
+double
+x_delta(double x)
+{
+	/* 
+	 * Calculate delta function value.
+	 */
+	return (x == 0) ? INFINITY : 0;
+}
+
+double
+x_nandelta(double x)
+{
+	/* 
+	 * Calculate modified delta function value.
+	 */
+	return (x == 0) ? NAN : 0;
+}
+
diff --git a/contrib/MathEval/xmath.h b/contrib/MathEval/xmath.h
index f086ad916a8fec8b85c4ffc3fb0ea008a6a2aed9..62003a22973a615acb295c8652ed855249ea7229 100644
--- a/contrib/MathEval/xmath.h
+++ b/contrib/MathEval/xmath.h
@@ -31,5 +31,9 @@ double x_acosh(double x);
 double x_atanh(double x);
 double x_actanh(double x);
 double x_rand(double x);
+double x_step(double x);
+double x_heavs(double x);
+double x_delta(double x);
+double x_nandelta(double x);
 
 #endif