diff --git a/bazaar/STEM4U/Integral.h b/bazaar/STEM4U/Integral.h
index 3dc33822c..c08f6eeeb 100644
--- a/bazaar/STEM4U/Integral.h
+++ b/bazaar/STEM4U/Integral.h
@@ -8,24 +8,30 @@ enum IntegralType {TRAPEZOIDAL, SIMPSON_1_3, SIMPSON_3_8, HERMITE_3, HERMITE_5};
 template <class Range, class T>
 T Integral(const Range &y, const Range &x, IntegralType type = TRAPEZOIDAL) {
 	ASSERT(x.size() == y.size());
-	ASSERT(x.size() > 1);
+	if (y.size() <= 1)
+		return 0;
 	
 	T ret = 0;
 	size_t n = x.size();
+	if(type == SIMPSON_1_3 && n == 2)
+		type = TRAPEZOIDAL;
+	else if (type == SIMPSON_3_8) {
+		if (n == 2)
+			type = TRAPEZOIDAL;
+		else if (n == 3)
+			type = SIMPSON_1_3;
+	}
+	
 	if (type == TRAPEZOIDAL) {
 		for (int i = 1; i < n; ++i)
 			ret += Avg(y(i), y(i-1))*(x(i) - x(i-1));
 	} else if (type == SIMPSON_1_3) {
-		if (n < 3)
-			return Null;
 		int i;
 		for (i = 2; i < n; i += 2)
 			ret += (x[i] - x[i-2])/6.*(y[i-2] + 4*y[i-1] + y[i]);
 		if (i == n)
 			ret += Avg(y(n-1), y(n-2))*(x(n-1) - x(n-2));
 	} else if (type == SIMPSON_3_8) {
-		if (n < 4)
-			return Null;
 		int i;
 		for (i = 3; i < n; i += 3)
 			ret += (x[i] - x[i-3])/8.*(y[i-3] + 3*y[i-2] + 3*y[i-1] + y[i]);
@@ -61,12 +67,25 @@ T Integral(Range &y, T dx, IntegralType type = TRAPEZOIDAL) {
 		return 0;
 	
 	Eigen::Index n = y.size();
+	if(type == SIMPSON_1_3 && n == 2)
+		type = TRAPEZOIDAL;
+	else if (type == SIMPSON_3_8) {
+		if (n == 2)
+			type = TRAPEZOIDAL;
+		else if (n == 3)
+			type = SIMPSON_1_3;
+	} else if (type == HERMITE_3 && n == 2) 
+		type = TRAPEZOIDAL;
+	else if (type == HERMITE_5) {
+		if (n == 2)
+			type = TRAPEZOIDAL;
+		else if (n == 3)
+			type = HERMITE_3;
+	}	
 	
 	if (type == TRAPEZOIDAL) 
 		return (y.segment(1, n-2).sum() + (y(0) + y(n-1))/2)*dx;
 	else if (type == SIMPSON_1_3) {
-		if (n < 3)
-			return Null;
 		T ret0 = 0;
 		if ((n-1)%2) {
 			ret0 = Avg(y(n-1), y(n-2))*dx;
@@ -74,8 +93,6 @@ T Integral(Range &y, T dx, IntegralType type = TRAPEZOIDAL) {
 		}
 		return ret0 + Calc1_3(y, dx, n);
 	} else if (type == SIMPSON_3_8) {
-		if (n < 4)
-			return Null;
 		T ret0 = 0;
 		int rem = (n-1)%3;
 		if (rem == 2) {
@@ -94,13 +111,11 @@ T Integral(Range &y, T dx, IntegralType type = TRAPEZOIDAL) {
 		}
 		return ret0 + ret*dx*3./8;
 	} else if (type == HERMITE_3) {
-		if (n < 3)
-			return Null;
 		return (y.segment(1, n-2).sum() + (y(0) + y(n-1))/2)*dx
 			 + dx/24.*(3*y[0] - 4*y[1] + y[2] + y[n-3] - 4*y[n-2] + 3*y[n-1]);
 	} else if (type == HERMITE_5) {
-		if (n < 5)
-			return Null;
+		if (n == 4)
+			return ((y[0] + y[3])/2. + y[1] + y[2])*dx;
 		return (y.segment(1, n-2).sum() + (y(0) + y(n-1))/2)*dx
 			 - dx/144.*(25*(y[0] - y[n-1]) - 48*(y[1] - y[n-2]) + 36*(y[2] - y[n-3])
 			 		  - 16*(y[3] - y[n-4]) +  3*(y[4] - y[n-5]));