diff --git a/romspline/__init__.py b/romspline/__init__.py
index 3c06c0e0667548778621059c2fc5355c0c16d48a..74bc810deb8c2710b3b5ca756edfe01b7728dbae 100644
--- a/romspline/__init__.py
+++ b/romspline/__init__.py
@@ -69,7 +69,7 @@ class _ImportStates(object):
       import numpy as np
     except:
       # Raise an exception if numpy can't be imported
-      raise Exception, "Error: Cannot import `NumPy` module."
+      raise Exception("Error: Cannot import `NumPy` module.")
     
     # Try importing scipy.interpolate.UnivariateSpline class
     try:
@@ -77,14 +77,14 @@ class _ImportStates(object):
     except:
       # Raise an exception if UnivariateSpline can't be imported
       # This class is crucial to RomSpline
-      raise Exception, "Error: Cannot import `scipy.interpolate.UnivariateSpline` class."
+      raise Exception("Error: Cannot import `scipy.interpolate.UnivariateSpline` class.")
       
     # Try importing h5py module
     try:
       import h5py
       self._H5PY = True
     except:
-      print "Warning: Cannot import `h5py` module. File I/O features will be limited to text formats."
+      print("Warning: Cannot import `h5py` module. File I/O features will be limited to text formats.")
       self._H5PY = False
     
     # Try importing matplotlib module
@@ -92,7 +92,7 @@ class _ImportStates(object):
       import matplotlib.pyplot as plt
       self._MATPLOTLIB = True
     except:
-      print "Warning: Cannot import `matplotlib.pyplot` module."
+      print("Warning: Cannot import `matplotlib.pyplot` module.")
       self._MATPLOTLIB = False
     
     # Try importing futures module
@@ -100,7 +100,7 @@ class _ImportStates(object):
       from concurrent.futures import ProcessPoolExecutor, wait, as_completed
       self._FUTURES = True
     except:
-      print "Warning: Cannot import `futures` module. Try running `pip install futures` to install."
+      print("Warning: Cannot import `futures` module. Try running `pip install futures` to install.")
       self._FUTURES = False
     
     # Try importing multiprocessing module
@@ -108,7 +108,7 @@ class _ImportStates(object):
       from multiprocessing import cpu_count
       self._MP = True
     except:
-      print "Warning: Cannot import `multiprocessing` module."
+      print("Warning: Cannot import `multiprocessing` module.")
       self._MP = False
     
     # If can import both futures and multiprocessing modules 
@@ -118,7 +118,7 @@ class _ImportStates(object):
     if self._FUTURES and self._MP:
       self._PARALLEL = True
     else:
-      print "Warning: Parallel computation options will be unavailable."
+      print("Warning: Parallel computation options will be unavailable.")
       self._PARALLEL = False
 
 state = _ImportStates()
@@ -128,12 +128,12 @@ state = _ImportStates()
 # Import submodules #
 #####################
 
-from greedy import *            # For building reduced-order splines
-from convergence import *       # For studying convergence
-from random_seeds import *      # For studying the effect of seed points on reduced data sizes
-from cross_validation import *  # For estimating (global) interpolation errors
-from build_spline import *      # Convenience module for bulding reduced-order spline
+from .greedy import *            # For building reduced-order splines
+from .convergence import *       # For studying convergence
+from .random_seeds import *      # For studying the effect of seed points on reduced data sizes
+from .cross_validation import *  # For estimating (global) interpolation errors
+from .build_spline import *      # Convenience module for bulding reduced-order spline
                                 # with a global interpolation error estimate from cross-validation
-from example import *         # Built-in function for testing and demonstration purposes
-from regression import *        # Regression testing
+from .example import *         # Built-in function for testing and demonstration purposes
+from .regression import *        # Regression testing
 
diff --git a/romspline/build_spline.py b/romspline/build_spline.py
index efb2836a4806025f20e807a6c53380cf70c8a120..54f95d04ead296ab8992aad3cbaff7cd65a69793 100644
--- a/romspline/build_spline.py
+++ b/romspline/build_spline.py
@@ -1,10 +1,10 @@
-from __init__ import state
+from .__init__ import state
 
 if state._MATPLOTLIB:
   import matplotlib.pyplot as plt
 
 import numpy as np
-import greedy, random_seeds, cross_validation
+from . import greedy, random_seeds, cross_validation
 
 
 
@@ -77,7 +77,7 @@ def build_spline(x, y, tol=1e-6, deg=None, rel=False, seeds=None, small=True, cv
   # Build a reduced-order spline with given specifications (deg, tol, etc.)
   if small is False:
     if verbose:
-      print "\nBuilding the reduced-order spline...",
+      print("\nBuilding the reduced-order spline...", end=' ')
     if deg is None:
       deg = 5
     spline = greedy.ReducedOrderSpline(x, y, tol=tol, deg=deg, rel=rel, seeds=seeds, verbose=verbose)
@@ -85,7 +85,7 @@ def build_spline(x, y, tol=1e-6, deg=None, rel=False, seeds=None, small=True, cv
   
   # Otherwise, sample the seed values to find a small reduced-order spline
   else:
-    print "\nFinding a smallest reduced-order data set..."
+    print("\nFinding a smallest reduced-order data set...")
     if small is True:  # Default number of seed sets to sample is 10
       num = 10
     else:              # Otherwise, `small` is the number of sets of seed points to sample
@@ -96,7 +96,7 @@ def build_spline(x, y, tol=1e-6, deg=None, rel=False, seeds=None, small=True, cv
     
   
   if cv is not False:
-    print "\nPerforming Monte Carlo K-fold cross-validation..."
+    print("\nPerforming Monte Carlo K-fold cross-validation...")
     if cv is True:  # Default number of Monte Carlo K-fold cross-validation studies to perform
       num = 10
     else:
diff --git a/romspline/convergence.py b/romspline/convergence.py
index 752c8a5a8489b61502d815bf16665be6a7105954..9b9f8aa02db3973cad9411cbfc3d1c494d52e22a 100644
--- a/romspline/convergence.py
+++ b/romspline/convergence.py
@@ -1,11 +1,11 @@
-from __init__ import state
+from .__init__ import state
 
 if state._MATPLOTLIB:
   import matplotlib.pyplot as plt
 
 import numpy as np
-import greedy
-from helpers import Linfty
+from . import greedy
+from .helpers import Linfty
 
 
 
@@ -203,7 +203,7 @@ class Convergence(object):
           return ax
       
     else:
-      print "No data to plot. Run `make` method."
+      print("No data to plot. Run `make` method.")
   
   
   def _plot_Linfty_errors(self, ax):
@@ -234,7 +234,7 @@ class Convergence(object):
           return ax
         
     else:
-      print "No data to plot. Run `make` method."
+      print("No data to plot. Run `make` method.")
   
   
   def _plot_sizes(self, ax, axes=None):
@@ -308,6 +308,6 @@ class Convergence(object):
           return ax
         
     else:
-      print "No data to plot. Run `make` method."
+      print("No data to plot. Run `make` method.")
 
 
diff --git a/romspline/cross_validation.py b/romspline/cross_validation.py
index 533e5d81eeb12afe3d2446cf4cb21696054e05b5..a702387df64808bbdf6bd82ee88dc46f91358f06 100644
--- a/romspline/cross_validation.py
+++ b/romspline/cross_validation.py
@@ -1,13 +1,13 @@
-from __init__ import state
+from .__init__ import state
 
 if state._MATPLOTLIB:
   import matplotlib.pyplot as plt
 if state._PARALLEL:
-  from __init__ import ProcessPoolExecutor, wait, as_completed, cpu_count
+  from .__init__ import ProcessPoolExecutor, wait, as_completed, cpu_count
 
 import numpy as np
-import greedy
-from helpers import *
+from . import greedy
+from .helpers import *
 
 
 
@@ -226,7 +226,7 @@ class CrossValidation(object):
     for nn in range(n):
       
       if verbose and not (nn+1)%10:
-        print "Trials completed:", nn+1
+        print("Trials completed:", nn+1)
       
       self.Kfold(x, y, K=K, parallel=parallel, random=random)
       self.all_args[nn] = self.args
@@ -327,7 +327,7 @@ class CrossValidation(object):
       if ax is None:
         fig, ax = plt.subplots(nrows=1, ncols=1)
       
-      x_plot = range(len(self._partitions))
+      x_plot = list(range(len(self._partitions)))
       ax.plot(x_plot, self.errors, color=color, marker=marker, linestyle=linestyle)
       
       ax.set_xlabel('Partition')
@@ -343,7 +343,7 @@ class CrossValidation(object):
           return ax
       
     else:
-      print "No data attributes to plot."
+      print("No data attributes to plot.")
   
   
   def plot_monte_carlo_errors(self, x=None, n=20, axes='plot', ax=None, show=True, color='k', marker='.'):
@@ -409,6 +409,6 @@ class CrossValidation(object):
           return ax
       
     else:
-      print "No data attributes to plot."
+      print("No data attributes to plot.")
 
 
diff --git a/romspline/example.py b/romspline/example.py
new file mode 100644
index 0000000000000000000000000000000000000000..8c98474a54804e11a4da48d0152cb5cdecffd839
--- /dev/null
+++ b/romspline/example.py
@@ -0,0 +1,134 @@
+import numpy as np
+
+
+################################
+#   Class for generating some  #
+#   test data for showing how  #
+#    to use the code in the    #
+#   Jupyter/IPython notebooks  #
+################################
+
+class TestData(object):
+  """Generate the test data used as in example IPython notebooks 
+    for demonstrating the construction, properties, and errors 
+    of a reduced-order spline interpolant.
+  """
+  
+  def __init__(self, num=4001, noise=0., uv=0.):
+    """Create a TestData object.
+    
+    Input
+    -----
+      num   -- number of samples to evaluate the function
+               in domain [-1,1]
+      noise -- amplitude of stochastic fluctuations added
+               to smooth function values
+               (default is 0.)
+      uv    -- amplitude of high-frequency (i.e., ultra-violet)
+               features added to smooth function values
+               (default is 0.)
+    
+    Attributes
+    ----------
+      x -- samples
+      y -- values of sampled function
+    
+    """
+    
+    # Generate test data
+    self.x = np.linspace(-1, 1, num)
+    self.y = self.f(self.x, noise=noise, uv=uv)
+  
+  
+  def f(self, x, noise=0., uv=0.):
+    """Function to sample for reduced-order spline examples
+    
+    Inputs
+    ------
+      x     -- values to sample the (smooth) function
+      noise -- amplitude of stochastic fluctuations added
+               to smooth function values
+               (default is 0.)
+      uv    -- amplitude of high-frequency (i.e., ultra-violet)
+               features added to smooth function values
+               (default is 0.)
+    
+    Output
+    ------
+      sampled function values
+    
+    Comments
+    --------
+    The function being evaluated is
+    
+      f(x) = 100.*( (1.+x) * sin(5.*(x-0.2)**2) 
+                    + exp(-(x-0.5)**2/2./0.01) * sin(100*x) 
+                  )
+    """
+    
+    # Validate inputs
+    x = np.asarray(x)
+    
+    # Return smooth function values
+    ans = 100.*( (x+1.)*np.sin(5.*(x-0.2)**2) + np.exp(-(x-0.5)**2/2./0.01)*np.sin(100*x) )
+    
+    # Return smooth function values with high-frequency (UV) features
+    if uv != 0.:
+      assert type(uv) in [float, int], "Expecting integer or float type."
+      ans += float(uv)*self.uv(x)
+    
+    # Return smooth function values with stochastic noise
+    if noise != 0.:
+      assert type(noise) in [float, int], "Expecting integer or float type."
+      ans += float(noise)*np.random.randn(len(x))
+    
+    return ans
+  
+  
+  def dfdx(self, x):
+    """Analytic derivative of f
+    
+    Inputs
+    ------
+      x -- values to sample the derivative of 
+           the function f(x) (see self.f method)
+    
+    Outputs
+    -------
+      ans -- values of analytically calculated 
+             derivative of the function
+    
+    """
+    x = np.asarray(x)
+    
+    a = 10.*(-0.2+x)*(1.+x)*np.cos(5.*(-0.2 + x)**2)
+    b = 100.*np.exp(-50.*(-0.5+x)**2)*np.cos(100.*x)
+    c = np.sin(5.*(-0.2+x)**2)
+    d = -100.*np.exp(-50.*(-0.5+x)**2)*(-0.5+x)*np.sin(100.*x)
+    
+    ans = 100.*(a+b+c+d)
+    
+    return ans
+  
+  
+  def uv(self, x, width=20):
+    """Generate high-frequency oscillations
+    
+    Inputs
+    ------
+      x     -- values to sample the high-frequency
+               oscillations
+      width -- number of samples corresponding to
+               the period of the high-frequency
+               oscillations
+               (default is 20)
+    
+    Outputs
+    -------
+      array of high-frequency oscillating values
+    
+    """
+    X = x[width] - x[0]
+    return np.sin(len(x)/X * x)
+
+
diff --git a/romspline/greedy.py b/romspline/greedy.py
index b7721ac5b771c335d8dc7a54722183d9525e717b..8cb0285271877659518ed7061a2a5e6da2322b33 100644
--- a/romspline/greedy.py
+++ b/romspline/greedy.py
@@ -1,4 +1,4 @@
-from __init__ import _ImportStates
+from .__init__ import _ImportStates
 
 state = _ImportStates()
 if state._MATPLOTLIB:
@@ -189,10 +189,10 @@ class ReducedOrderSpline(object):
       if abs:
         errors = np.abs(errors)
     else:
-      raise Exception, "No spline interpolant to compare against. Run the greedy method."
+      raise Exception("No spline interpolant to compare against. Run the greedy method.")
     
     if verbose:
-      print "Requested tolerance of {} met: ".format(self._tol), np.all(np.abs(errors) <= self.tol)
+      print("Requested tolerance of {} met: ".format(self._tol), np.all(np.abs(errors) <= self.tol))
     
     return errors
   
@@ -251,10 +251,10 @@ class ReducedOrderSpline(object):
             return ax
         
       else:
-        print "No data to plot. Run `greedy` method."
+        print("No data to plot. Run `greedy` method.")
     
     else:
-      print "Error: matplotlib.pyplot module not imported."
+      print("Error: matplotlib.pyplot module not imported.")
   
   
   def write(self, file, slim=False):
@@ -289,7 +289,7 @@ class ReducedOrderSpline(object):
     """
     
     if not state._H5PY:
-      print "h5py module not imported. Try writing data to text (.txt) format."
+      print("h5py module not imported. Try writing data to text (.txt) format.")
       return
     
     # If file is an HDF5 file or group descriptor...
@@ -310,12 +310,12 @@ class ReducedOrderSpline(object):
             fp = h5py.File(file, 'w')
             isopen = True
           except:
-            raise Exception, "Could not open file for writing."
+            raise Exception("Could not open file for writing.")
           if isopen:
             self._write(fp, slim=slim)
             fp.close()
         else:
-          print "Error: h5py module is not imported. Try writing data to text (.txt) format."
+          print("Error: h5py module is not imported. Try writing data to text (.txt) format.")
           return
       
       # Text format
@@ -364,7 +364,7 @@ class ReducedOrderSpline(object):
       if not slim:
         descriptor.create_dataset('errors', data=self.errors, dtype='double', compression='gzip', shuffle=True)
     else:
-      raise Exception, "Descriptor not recognized."
+      raise Exception("Descriptor not recognized.")
   
   
   def read(self, file, group=None):
@@ -418,21 +418,21 @@ def readSpline(file, group=None):
         fp = h5py.File(file, 'r')
         isopen = True
       except:
-        raise Exception, "Could not open file for reading."
+        raise Exception("Could not open file for reading.")
       if isopen:
         gp = fp[group] if group else fp
         deg = gp['deg'][()]
         tol = gp['tol'][()]
         X = gp['X'][:]
         Y = gp['Y'][:]
-        if hasattr(gp, 'errors') or 'errors' in gp.keys():
+        if hasattr(gp, 'errors') or 'errors' in list(gp.keys()):
           errors = gp['errors'][:]
         else:
           errors = []
         fp.close()
         _made = True
     else:
-      print "Error: h5py module is not imported."
+      print("Error: h5py module is not imported.")
       return
   
   # Text format
@@ -449,7 +449,7 @@ def readSpline(file, group=None):
         errs_isopen = False
       isopen = True
     except:
-      raise IOError, "Could not open file(s) for reading."
+      raise IOError("Could not open file(s) for reading.")
     
     if isopen:
       deg = int(fp_deg.read())
@@ -492,7 +492,7 @@ def readSpline(file, group=None):
     spline._made = _made
     return spline
   else:
-    raise Exception, "Reduced-order spline interpolant could not be constructed from file."
+    raise Exception("Reduced-order spline interpolant could not be constructed from file.")
 
 
 
@@ -537,8 +537,8 @@ def _seed(x, deg=5, seeds=None):
 def _greedy(x, y, tol=1e-6, rel=False, deg=5, verbose=False, seeds=None):
   """Greedy algorithm for building a reduced-order spline"""
   if verbose:
-    print "\nSize", "\t", "Error"
-    print "="*13
+    print("\nSize", "\t", "Error")
+    print("="*13)
   
   if rel:
     ymax = np.max(np.abs(y))
@@ -570,7 +570,7 @@ def _greedy(x, y, tol=1e-6, rel=False, deg=5, verbose=False, seeds=None):
     
     # Print to screen, if requested
     if verbose:
-      print ctr, "\t", errors[-1]
+      print(ctr, "\t", errors[-1])
     
     # Check if greedy error is below tolerance and exit if so
     if errors[-1] < _tol:
diff --git a/romspline/h5splinetotxt.py b/romspline/h5splinetotxt.py
index d9c3d5aa9a2020413a3e262b5194515e72564014..44a67e9bdba6ecac96d16649407a15f805849edd 100644
--- a/romspline/h5splinetotxt.py
+++ b/romspline/h5splinetotxt.py
@@ -17,4 +17,4 @@ coords = f[coords_dataset]
 vals = s_h5(coords)
 
 for t,f in zip(coords,vals):
-    print "%.19g\t%.19g" % (t,f)
+    print("%.19g\t%.19g" % (t,f))
diff --git a/romspline/helpers.py b/romspline/helpers.py
index 6510ccb27bb654be699658ac60da223a85a96d81..4abce4a80ee55e80b9bafafdb15d1c21e9aa4a59 100644
--- a/romspline/helpers.py
+++ b/romspline/helpers.py
@@ -1,4 +1,4 @@
-from __init__ import state
+from .__init__ import state
 import numpy as np
 
 
@@ -43,7 +43,7 @@ def random_partitions(n, K):
   assert n >= K, "Number of partitions must not exceed number of samples."
   
   # Make array with unique random integers
-  rand = np.random.choice(range(n), n, replace=False)
+  rand = np.random.choice(list(range(n)), n, replace=False)
   
   # Split into K (nearly) equal partitions
   return [np.sort(rand[pp]) for pp in partitions(n, K)]
diff --git a/romspline/random_seeds.py b/romspline/random_seeds.py
index da212ef748b7766d8e3874ba7576fc82441868bd..9a79679a5bb9359210ab451efc5a0ea2acb8a1bb 100644
--- a/romspline/random_seeds.py
+++ b/romspline/random_seeds.py
@@ -1,4 +1,4 @@
-from __init__ import state
+from .__init__ import state
 
 if state._MATPLOTLIB:
   import matplotlib.pyplot as plt
@@ -7,7 +7,7 @@ if state._PARALLEL:
   from multiprocessing import cpu_count
 
 import numpy as np
-import greedy
+from . import greedy
 
 
 
@@ -139,7 +139,7 @@ class RandomSeeds(object):
     
     if (parallel is False) or (state._PARALLEL is False):
       for nn in range(Nseeds):
-        seeds = np.sort( np.random.choice(range(len(x)), self._deg+1, replace=False) )
+        seeds = np.sort( np.random.choice(list(range(len(x))), self._deg+1, replace=False) )
         self.errors[nn], self.sizes[nn], self.seeds[nn], Xs = _randomSeeds(x, y, seeds, tol=self._tol, rel=self._rel, deg=self._deg)
         self.Xs.append(Xs)
     elif state._PARALLEL is True:
@@ -158,7 +158,7 @@ class RandomSeeds(object):
       # Build reduced-order splines for many random sets of seed points
       output = []
       for nn in range(Nseeds):
-        seeds = np.sort( np.random.choice(range(len(x)), self._deg+1, replace=False) )
+        seeds = np.sort( np.random.choice(list(range(len(x))), self._deg+1, replace=False) )
         output.append(executor.submit(_randomSeeds, x, y, seeds, tol=self._tol, rel=self._rel, deg=self._deg))
       
       # Gather the results as they complete
@@ -219,7 +219,7 @@ class RandomSeeds(object):
           return ax
     
     else:
-      print "No data to plot. Run make method or set _made attribute to True."
+      print("No data to plot. Run make method or set _made attribute to True.")
     
   
   def plot_samples(self, x, y, ax=None, show=True, xlabel='$x$', ylabel='Size of reduced data, $n$'):
@@ -228,7 +228,7 @@ class RandomSeeds(object):
         fig, ax = plt.subplots(nrows=1, ncols=1)
       
       # Assemble all reduced samples into a 2d array
-      longest = max(map(len, self.Xs))
+      longest = max(list(map(len, self.Xs)))
       Xs = max(x)*np.ones((len(self.Xs),longest), dtype='double')
       for ii, kk in enumerate(self.Xs):
         Xs[ii][:len(kk)] = kk
@@ -256,7 +256,7 @@ class RandomSeeds(object):
           return ax
     
     else:
-      print "No data to plot. Run `make` method or set _made attribute to True."
+      print("No data to plot. Run `make` method or set _made attribute to True.")
 
 
 
@@ -294,14 +294,14 @@ def small_spline(x, y, num, tol=1e-6, deg=None, rel=False, parallel=True, verbos
   """
   
   if deg is None:                  # Use all integers in [1,5]
-    degs = range(1, 6)
+    degs = list(range(1, 6))
   else:
     if len(np.shape(deg)) == 1:    # deg is a list or array
       degs = deg
     elif len(np.shape(deg)) == 0:  # deg is a number
       degs = [deg]
     else:
-      raise Exception, "Input for `deg` option not recognized."
+      raise Exception("Input for `deg` option not recognized.")
   for dd in degs:
     assert (dd in range(1,6)), "Expecting degree(s) to be one or more integers in [1,5]."
   
@@ -309,7 +309,7 @@ def small_spline(x, y, num, tol=1e-6, deg=None, rel=False, parallel=True, verbos
   
   for ii, dd in enumerate(degs):
     if verbose:
-      print "Smallest spline for degree {} is...".format(dd),
+      print("Smallest spline for degree {} is...".format(dd), end=' ')
     
     # Sample from the set of all possible initial 
     # seeds for this polynomial degree
@@ -317,7 +317,7 @@ def small_spline(x, y, num, tol=1e-6, deg=None, rel=False, parallel=True, verbos
     rand.make(x, y, num, parallel=parallel)
     
     if verbose:  # Print smallest size found in sample
-      print int(np.min(rand.sizes))
+      print(int(np.min(rand.sizes)))
     
     # Find smallest spline in the sample for this polynomial degree
     imin = np.argmin(rand.sizes)
diff --git a/romspline/regression.py b/romspline/regression.py
index a88dcb531d93304d7038f16d36992bb95133959f..9e3f6dcfb5840e5fb985f98596e983e980f63510 100644
--- a/romspline/regression.py
+++ b/romspline/regression.py
@@ -3,13 +3,13 @@ Test the reduced-order spline greedy algorithm outputs with
 previously generated data to make sure the code still works correctly.
 """
 
-from __init__ import state
+from .__init__ import state
 
 if state._H5PY:
   import h5py
 
 import numpy as np
-import greedy, example
+from . import greedy, example
 
 
 def regression():
@@ -36,29 +36,29 @@ def regression():
   
   # Build reduced-order spline interpolant for the
   # test data using options found in regressionData.h5.
-  print "Building reduced-order spline...",
+  print("Building reduced-order spline...", end=' ')
   test_spline = greedy.ReducedOrderSpline(test.x, test.y, deg=deg, tol=tol)
-  print "Done"
+  print("Done")
   
   # Perform a few checks on the greedy algorithm
-  print "Testing greedy algorithm:"
+  print("Testing greedy algorithm:")
   
-  print "  Comparing reduced data size...",
+  print("  Comparing reduced data size...", end=' ')
   if test_spline.size == len(X):
-    print "Passed"
+    print("Passed")
   else:
-    print "Failed [Size of reduced data ({}) does not equal {}]".format(test_spline.size, len(X))
+    print("Failed [Size of reduced data ({}) does not equal {}]".format(test_spline.size, len(X)))
   
-  print "  Comparing selected array indices...",
+  print("  Comparing selected array indices...", end=' ')
   if np.all([test_spline.args[ii] == args[ii] for ii in range(len(args))]):
-    print "Passed"
+    print("Passed")
   else:
-    print "Failed"
+    print("Failed")
   
-  print "  Comparing greedy errors...",
+  print("  Comparing greedy errors...", end=' ')
   if np.all([np.isclose(test_spline.errors[ii], errors[ii]) for ii in range(len(errors))]):
-    print "Passed"
+    print("Passed")
   else:
-    print "Failed"
+    print("Failed")