diff --git a/THANKS.txt b/THANKS.txt
index 0d277907b2d5..cee34775fd52 100644
--- a/THANKS.txt
+++ b/THANKS.txt
@@ -146,6 +146,7 @@ Johannes Ballé for the generalized normal distribution.
 Irvin Probst (ENSTA Bretagne) for pole placement.
 Ian Henriksen for Cython wrappers for BLAS and LAPACK
 Fukumu Tsutsumi for bug fixes.
+J.J. Green for interpolation bug fixes.
 
 
 Institutions
diff --git a/scipy/interpolate/rbf.py b/scipy/interpolate/rbf.py
index 45b876ad680a..53ab082c6b2a 100644
--- a/scipy/interpolate/rbf.py
+++ b/scipy/interpolate/rbf.py
@@ -45,12 +45,11 @@
 from __future__ import division, print_function, absolute_import
 
 import sys
+import numpy as np
 
-from numpy import (sqrt, log, asarray, newaxis, all, dot, exp, eye,
-                   float_)
 from scipy import linalg
-from scipy._lib.six import callable, get_method_function, \
-     get_function_code
+from scipy._lib.six import callable, get_method_function, get_function_code
+from scipy.special import xlogy
 
 __all__ = ['Rbf']
 
@@ -111,16 +110,16 @@ def euclidean_norm(x1, x2):
     """
 
     def _euclidean_norm(self, x1, x2):
-        return sqrt(((x1 - x2)**2).sum(axis=0))
+        return np.sqrt(((x1 - x2)**2).sum(axis=0))
 
     def _h_multiquadric(self, r):
-        return sqrt((1.0/self.epsilon*r)**2 + 1)
+        return np.sqrt((1.0/self.epsilon*r)**2 + 1)
 
     def _h_inverse_multiquadric(self, r):
-        return 1.0/sqrt((1.0/self.epsilon*r)**2 + 1)
+        return 1.0/np.sqrt((1.0/self.epsilon*r)**2 + 1)
 
     def _h_gaussian(self, r):
-        return exp(-(1.0/self.epsilon*r)**2)
+        return np.exp(-(1.0/self.epsilon*r)**2)
 
     def _h_linear(self, r):
         return r
@@ -132,9 +131,7 @@ def _h_quintic(self, r):
         return r**5
 
     def _h_thin_plate(self, r):
-        result = r**2 * log(r)
-        result[r == 0] = 0  # the spline is zero at zero
-        return result
+        return xlogy(r**2, r)
 
     # Setup self._function and do smoke test on initial r
     def _init_function(self, r):
@@ -186,10 +183,10 @@ def _init_function(self, r):
         return a0
 
     def __init__(self, *args, **kwargs):
-        self.xi = asarray([asarray(a, dtype=float_).flatten()
+        self.xi = np.asarray([np.asarray(a, dtype=np.float_).flatten()
                            for a in args[:-1]])
         self.N = self.xi.shape[-1]
-        self.di = asarray(args[-1]).flatten()
+        self.di = np.asarray(args[-1]).flatten()
 
         if not all([x.size == self.di.size for x in self.xi]):
             raise ValueError("All arrays must be equal length.")
@@ -198,7 +195,11 @@ def __init__(self, *args, **kwargs):
         r = self._call_norm(self.xi, self.xi)
         self.epsilon = kwargs.pop('epsilon', None)
         if self.epsilon is None:
-            self.epsilon = r.mean()
+            # default epsilon is the "the average distance between nodes"
+            dim = self.xi.shape[0]
+            ximax = np.amax(self.xi, axis=1)
+            ximin = np.amin(self.xi, axis=1)
+            self.epsilon = np.power(np.prod(ximax - ximin)/self.N, 1.0/dim)
         self.smooth = kwargs.pop('smooth', 0.0)
 
         self.function = kwargs.pop('function', 'multiquadric')
@@ -209,23 +210,23 @@ def __init__(self, *args, **kwargs):
         for item, value in kwargs.items():
             setattr(self, item, value)
 
-        self.A = self._init_function(r) - eye(self.N)*self.smooth
+        self.A = self._init_function(r) - np.eye(self.N)*self.smooth
         self.nodes = linalg.solve(self.A, self.di)
 
     def _call_norm(self, x1, x2):
         if len(x1.shape) == 1:
-            x1 = x1[newaxis, :]
+            x1 = x1[np.newaxis, :]
         if len(x2.shape) == 1:
-            x2 = x2[newaxis, :]
-        x1 = x1[..., :, newaxis]
-        x2 = x2[..., newaxis, :]
+            x2 = x2[np.newaxis, :]
+        x1 = x1[..., :, np.newaxis]
+        x2 = x2[..., np.newaxis, :]
         return self.norm(x1, x2)
 
     def __call__(self, *args):
-        args = [asarray(x) for x in args]
+        args = [np.asarray(x) for x in args]
         if not all([x.shape == y.shape for x in args for y in args]):
             raise ValueError("Array lengths must be equal")
         shp = args[0].shape
-        xa = asarray([a.flatten() for a in args], dtype=float_)
+        xa = np.asarray([a.flatten() for a in args], dtype=np.float_)
         r = self._call_norm(xa, self.xi)
-        return dot(self._function(r), self.nodes).reshape(shp)
+        return np.dot(self._function(r), self.nodes).reshape(shp)
diff --git a/scipy/interpolate/tests/test_rbf.py b/scipy/interpolate/tests/test_rbf.py
index a20215919422..fda5c1f01174 100644
--- a/scipy/interpolate/tests/test_rbf.py
+++ b/scipy/interpolate/tests/test_rbf.py
@@ -15,48 +15,36 @@
 
 
 def check_rbf1d_interpolation(function):
-    """Check that the Rbf function interpolates through the nodes (1D)"""
-    olderr = np.seterr(all="ignore")
-    try:
-        x = linspace(0,10,9)
-        y = sin(x)
-        rbf = Rbf(x, y, function=function)
-        yi = rbf(x)
-        assert_array_almost_equal(y, yi)
-        assert_almost_equal(rbf(float(x[0])), y[0])
-    finally:
-        np.seterr(**olderr)
+    # Check that the Rbf function interpolates through the nodes (1D)
+    x = linspace(0,10,9)
+    y = sin(x)
+    rbf = Rbf(x, y, function=function)
+    yi = rbf(x)
+    assert_array_almost_equal(y, yi)
+    assert_almost_equal(rbf(float(x[0])), y[0])
 
 
 def check_rbf2d_interpolation(function):
-    """Check that the Rbf function interpolates through the nodes (2D)"""
-    olderr = np.seterr(all="ignore")
-    try:
-        x = random.rand(50,1)*4-2
-        y = random.rand(50,1)*4-2
-        z = x*exp(-x**2-1j*y**2)
-        rbf = Rbf(x, y, z, epsilon=2, function=function)
-        zi = rbf(x, y)
-        zi.shape = x.shape
-        assert_array_almost_equal(z, zi)
-    finally:
-        np.seterr(**olderr)
+    # Check that the Rbf function interpolates through the nodes (2D).
+    x = random.rand(50,1)*4-2
+    y = random.rand(50,1)*4-2
+    z = x*exp(-x**2-1j*y**2)
+    rbf = Rbf(x, y, z, epsilon=2, function=function)
+    zi = rbf(x, y)
+    zi.shape = x.shape
+    assert_array_almost_equal(z, zi)
 
 
 def check_rbf3d_interpolation(function):
-    """Check that the Rbf function interpolates through the nodes (3D)"""
-    olderr = np.seterr(all="ignore")
-    try:
-        x = random.rand(50,1)*4-2
-        y = random.rand(50,1)*4-2
-        z = random.rand(50,1)*4-2
-        d = x*exp(-x**2-y**2)
-        rbf = Rbf(x, y, z, d, epsilon=2, function=function)
-        di = rbf(x, y, z)
-        di.shape = x.shape
-        assert_array_almost_equal(di, d)
-    finally:
-        np.seterr(**olderr)
+    # Check that the Rbf function interpolates through the nodes (3D).
+    x = random.rand(50, 1)*4 - 2
+    y = random.rand(50, 1)*4 - 2
+    z = random.rand(50, 1)*4 - 2
+    d = x*exp(-x**2 - y**2)
+    rbf = Rbf(x, y, z, d, epsilon=2, function=function)
+    di = rbf(x, y, z)
+    di.shape = x.shape
+    assert_array_almost_equal(di, d)
 
 
 def test_rbf_interpolation():
@@ -67,31 +55,28 @@ def test_rbf_interpolation():
 
 
 def check_rbf1d_regularity(function, atol):
-    """Check that the Rbf function approximates a smooth function well away
-    from the nodes."""
-    olderr = np.seterr(all="ignore")
-    try:
-        x = linspace(0, 10, 9)
-        y = sin(x)
-        rbf = Rbf(x, y, function=function)
-        xi = linspace(0, 10, 100)
-        yi = rbf(xi)
-        #import matplotlib.pyplot as plt
-        #plt.figure()
-        #plt.plot(x, y, 'o', xi, sin(xi), ':', xi, yi, '-')
-        #plt.title(function)
-        #plt.show()
-        msg = "abs-diff: %f" % abs(yi - sin(xi)).max()
-        assert_(allclose(yi, sin(xi), atol=atol), msg)
-    finally:
-        np.seterr(**olderr)
+    # Check that the Rbf function approximates a smooth function well away
+    # from the nodes.
+    x = linspace(0, 10, 9)
+    y = sin(x)
+    rbf = Rbf(x, y, function=function)
+    xi = linspace(0, 10, 100)
+    yi = rbf(xi)
+    # import matplotlib.pyplot as plt
+    # plt.figure()
+    # plt.plot(x, y, 'o', xi, sin(xi), ':', xi, yi, '-')
+    # plt.plot(x, y, 'o', xi, yi-sin(xi), ':')
+    # plt.title(function)
+    # plt.show()
+    msg = "abs-diff: %f" % abs(yi - sin(xi)).max()
+    assert_(allclose(yi, sin(xi), atol=atol), msg)
 
 
 def test_rbf_regularity():
     tolerances = {
-        'multiquadric': 0.05,
-        'inverse multiquadric': 0.02,
-        'gaussian': 0.01,
+        'multiquadric': 0.1,
+        'inverse multiquadric': 0.15,
+        'gaussian': 0.15,
         'cubic': 0.15,
         'quintic': 0.1,
         'thin-plate': 0.1,
@@ -101,9 +86,30 @@ def test_rbf_regularity():
         yield check_rbf1d_regularity, function, tolerances.get(function, 1e-2)
 
 
+def check_rbf1d_stability(function):
+    # Check that the Rbf function with default epsilon is not subject 
+    # to overshoot.  Regression for issue #4523.
+    #
+    # Generate some data (fixed random seed hence deterministic) 
+    np.random.seed(1234)
+    x = np.linspace(0, 10, 50)
+    z = x + 4.0 * np.random.randn(len(x))
+
+    rbf = Rbf(x, z, function=function)
+    xi = np.linspace(0, 10, 1000)
+    yi = rbf(xi)
+
+    # subtract the linear trend and make sure there no spikes
+    assert_(np.abs(yi-xi).max() / np.abs(z-x).max() < 1.1)
+
+def test_rbf_stability():
+    for function in FUNCTIONS:
+        yield check_rbf1d_stability, function
+
+
 def test_default_construction():
-    """Check that the Rbf class can be constructed with the default
-    multiquadric basis function. Regression test for ticket #1228."""
+    # Check that the Rbf class can be constructed with the default
+    # multiquadric basis function. Regression test for ticket #1228.
     x = linspace(0,10,9)
     y = sin(x)
     rbf = Rbf(x, y)
@@ -112,7 +118,7 @@ def test_default_construction():
 
 
 def test_function_is_callable():
-    """Check that the Rbf class can be constructed with function=callable."""
+    # Check that the Rbf class can be constructed with function=callable.
     x = linspace(0,10,9)
     y = sin(x)
     linfunc = lambda x:x
@@ -122,8 +128,8 @@ def test_function_is_callable():
 
 
 def test_two_arg_function_is_callable():
-    """Check that the Rbf class can be constructed with a two argument
-    function=callable."""
+    # Check that the Rbf class can be constructed with a two argument
+    # function=callable.
     def _func(self, r):
         return self.epsilon + r