diff --git a/interpolation/interp.c b/interpolation/interp.c
index bb75e10..275735f 100644
--- a/interpolation/interp.c
+++ b/interpolation/interp.c
@@ -142,7 +142,18 @@ gsl_interp_eval (const gsl_interp * interp,
                  gsl_interp_accel * a)
 {
   double y;
-  int status = interp->type->eval (interp->state, xa, ya, interp->size, x, a, &y);
+  int status;
+
+  if (x < interp->xmin)
+    {
+      GSL_ERROR_VAL("interpolation error", GSL_EDOM, ya[0]);
+    }
+  else if (x > interp->xmax)
+    {
+      GSL_ERROR_VAL("interpolation error", GSL_EDOM, ya[interp->size - 1]);
+    }
+
+  status = interp->type->eval (interp->state, xa, ya, interp->size, x, a, &y);
 
   DISCARD_STATUS(status);
 
@@ -176,7 +187,18 @@ gsl_interp_eval_deriv (const gsl_interp * interp,
                        gsl_interp_accel * a)
 {
   double dydx;
-  int status = interp->type->eval_deriv (interp->state, xa, ya, interp->size, x, a, &dydx);
+  int status;
+
+  if (x < interp->xmin)
+    {
+      GSL_ERROR_VAL("interpolation error", GSL_EDOM, 0.0);
+    }
+  else if (x > interp->xmax)
+    {
+      GSL_ERROR_VAL("interpolation error", GSL_EDOM, 0.0);
+    }
+
+  status = interp->type->eval_deriv (interp->state, xa, ya, interp->size, x, a, &dydx);
 
   DISCARD_STATUS(status);
 
@@ -210,7 +232,18 @@ gsl_interp_eval_deriv2 (const gsl_interp * interp,
                         gsl_interp_accel * a)
 {
   double d2;
-  int status = interp->type->eval_deriv2 (interp->state, xa, ya, interp->size, x, a, &d2);
+  int status;
+
+  if (x < interp->xmin)
+    {
+      GSL_ERROR_VAL("interpolation error", GSL_EDOM, 0.0);
+    }
+  else if (x > interp->xmax)
+    {
+      GSL_ERROR_VAL("interpolation error", GSL_EDOM, 0.0);
+    }
+
+  status = interp->type->eval_deriv2 (interp->state, xa, ya, interp->size, x, a, &d2);
 
   DISCARD_STATUS(status);
 
@@ -247,7 +280,18 @@ gsl_interp_eval_integ (const gsl_interp * interp,
                        gsl_interp_accel * acc)
 {
   double result;
-  int status = interp->type->eval_integ (interp->state, xa, ya, interp->size, acc, a, b, &result);
+  int status;
+
+  if (a > b || a < interp->xmin || b > interp->xmax)
+    {
+      GSL_ERROR_VAL("interpolation error", GSL_EDOM, 0.0);
+    }
+  else if(a == b)
+    {
+      return 0.0;
+    }
+
+  status = interp->type->eval_integ (interp->state, xa, ya, interp->size, acc, a, b, &result);
 
   DISCARD_STATUS(status);