diff --git a/odespy/PyDSTool.py b/odespy/PyDSTool.py
index 3db4e76..a2fb76c 100644
--- a/odespy/PyDSTool.py
+++ b/odespy/PyDSTool.py
@@ -1,7 +1,9 @@
 # Author: Liwei Wang
 """
 """
-from solvers import *
+from __future__ import print_function
+from __future__ import absolute_import
+from .solvers import *
 import numpy as np
 
 class Pyds(Solver):
@@ -24,10 +26,10 @@ class Pyds(Solver):
 
     def initialize(self):
         try:
-            import PyDSTool
+            from . import PyDSTool
         except ImportError:
-            raise ImportError,'''
-        PyDSTool is not installed - required for solvers from PyDSTool'''
+            raise ImportError('''
+        PyDSTool is not installed - required for solvers from PyDSTool''')
 
     def solve(self, time_points, terminate=None):
         # Common parts as superclass
@@ -45,7 +47,7 @@ def solve(self, time_points, terminate=None):
         # through Python dictionaries with string keys.
                 
         # Start setting for PyDSTool
-        import PyDSTool
+        from . import PyDSTool
         neq, f, u0 = self.neq, self.f, self.U0
         
         # Initialize variables as trajectories in PyDSTOOL
@@ -148,7 +150,7 @@ def initialize_for_solve(self):
     method = Vode_pyds(f)
     method.set_initial_condition([0.,1.])
     u,t = method.solve(np.linspace(0.,10.,50))
-    print u
+    print(u)
     import scitools.std as st
     st.plot(t,u[:,0])
-    print max(u[:,0]-np.sin(t))
+    print(max(u[:,0]-np.sin(t)))
diff --git a/odespy/RungeKutta.py b/odespy/RungeKutta.py
index 73f99a1..ac39e8d 100644
--- a/odespy/RungeKutta.py
+++ b/odespy/RungeKutta.py
@@ -1,4 +1,6 @@
-from solvers import Solver, Adaptive
+from __future__ import print_function
+from __future__ import absolute_import
+from .solvers import Solver, Adaptive
 import numpy as np
 
 def _calculate_order_1_level(coefficients):
@@ -134,7 +136,7 @@ def middle(x,y,z):    # Auxilary function
         k = np.zeros((k_len, self.neq), self.dtype)  # intern stages
 
         if self.verbose > 0:
-            print 'advance solution in [%s, %s], h=%g' % (t_n, t_next, h)
+            print('advance solution in [%s, %s], h=%g' % (t_n, t_next, h))
 
         # Loop until next time point is reached
         while (abs(t - t_n) < abs(t_next - t_n)):
@@ -150,7 +152,7 @@ def middle(x,y,z):    # Auxilary function
 
             self.info['rejected'] += 1  # reduced below if accepted
             if self.verbose > 0:
-                print '  u(t=%g)=%g: ' % (t+h, u_new),
+                print('  u(t=%g)=%g: ' % (t+h, u_new), end=' ')
 
             # local error between 2 levels
             error = h*np.abs(np.dot(factors_error, k))
@@ -171,18 +173,18 @@ def middle(x,y,z):    # Auxilary function
                 self.info['rejected'] -= 1
 
                 if self.verbose > 0:
-                    print 'accepted, ',
+                    print('accepted, ', end=' ')
             else:
                 if self.verbose > 0:
-                    print 'rejected, ',
+                    print('rejected, ', end=' ')
 
             if self.verbose > 0:
-                print 'err=%s, ' % str(error),
+                print('err=%s, ' % str(error), end=' ')
                 if hasattr(self, 'u_exact') and callable(self.u_exact):
-                    print 'exact-err=%s, ' % \
-                          (np.asarray(self.u_exact(t+h))-u_new),
+                    print('exact-err=%s, ' % \
+                          (np.asarray(self.u_exact(t+h))-u_new), end=' ')
                 if h <= self.min_step:
-                    print 'h=min_step!! ',
+                    print('h=min_step!! ', end=' ')
 
 
            # Replace 0 values by 1e-16 since we will divide by error
@@ -209,7 +211,7 @@ def middle(x,y,z):    # Auxilary function
             h = min(h, t_next - t_intermediate[-1])
 
             if self.verbose > 0:
-                print 'new h=%g' % h
+                print('new h=%g' % h)
 
             if h == 0:
                 break
@@ -367,16 +369,16 @@ def validate_data(self):
         # Check for dimension of user-defined butcher table.
         array_shape = self.butcher_tableau.shape
         if len(array_shape) is not 2:
-            raise ValueError,'''
-        Illegal input! Your input butcher_tableau should be a 2d-array!'''
+            raise ValueError('''
+        Illegal input! Your input butcher_tableau should be a 2d-array!''')
         else:
             m,n = array_shape
             if m not in (n, n + 1):
-                raise ValueError, '''\
+                raise ValueError('''\
         The dimension of 2d-array <method_yours_array> should be:
         1. Either (n, n), --> For 1-level RungeKutta methods
         2. Or (n+1, n),   --> For 2-levels RungeKutta methods
-        The shape of your input array is (%d, %d).''' % (m,n)
+        The shape of your input array is (%d, %d).''' % (m,n))
         self._butcher_tableau = self.butcher_tableau
 
         # Check for user-defined order,
@@ -393,19 +395,19 @@ def validate_data(self):
             if array_shape[0] == array_shape[1] + 1:
                 # 2-level RungeKutta methods
                 if type(self.method_order) is int:
-                    raise ValueError, error_2level
+                    raise ValueError(error_2level)
                 try:
                     order1, order2 = self.method_order
                     if abs(order1-order2) != 1 or \
                             order1 < 1 or order2 < 1:
-                        raise ValueError, error_2level
+                        raise ValueError(error_2level)
                 except:
-                    raise ValueError,error_2level
+                    raise ValueError(error_2level)
             else:
                 # 1-level RungeKutta methods
                 if type(self.method_order) is not int or \
                         self.method_order < 1:
-                    raise ValueError,error_1level
+                    raise ValueError(error_1level)
             self._method_order = self.method_order
 
         else:   # method_order is not specified
@@ -418,7 +420,7 @@ def validate_data(self):
         for i in range(1,array_shape[1] - 1):
             if not np.allclose(self.butcher_tableau[i][0],\
                                sum(self.butcher_tableau[i][1:])):
-                raise ValueError, '''
+                raise ValueError('''
         Inconsistent data in Butcher_Tableau!
         In each lines of stage-coefficients, first number should be
         equal to the sum of other numbers.
@@ -426,6 +428,6 @@ def validate_data(self):
             a[i][0] == a[i][1] + a[i][2] + ... + a[i][K - 1]
             where 1 <= i <= K - 1
         Your input for line %d is :%s
-        ''' % (i,str(self.butcher_tableau[i]))
+        ''' % (i,str(self.butcher_tableau[i])))
 
         return True
diff --git a/odespy/__init__.py b/odespy/__init__.py
index 9212728..6d94bf8 100644
--- a/odespy/__init__.py
+++ b/odespy/__init__.py
@@ -5,6 +5,7 @@
 supported.  A wide range of numerical methods for ODEs are offered:
 
 """
+from __future__ import absolute_import
 
 # Insert tutorial from ../doc/src/odespy/odespy.rst
 
@@ -1255,13 +1256,13 @@ def f(u, t):
 with :math:`f(u,t)` implemented in Fortran.
 '''
 
-from solvers import *
-from RungeKutta import *
-from rkc import *
-from rkf45 import *
-from odepack import *
-from radau5 import *
-import problems
+from .solvers import *
+from .RungeKutta import *
+from .rkc import *
+from .rkf45 import *
+from .odepack import *
+from .radau5 import *
+from . import problems
 
 # Update doc strings with common info
 class_, doc_str, classname = None, None, None
@@ -1283,7 +1284,7 @@ def f(u, t):
 
 # Do not pollute namespace
 del class_, doc_str, classname, classnames, toc, typeset_toc, \
-    table_of_parameters, name, obj, inspect
+    table_of_parameters, inspect
 
 if __name__ == '__main__':
     from os.path import join
diff --git a/odespy/demos/demo_Lsodi_1.py b/odespy/demos/demo_Lsodi_1.py
index 84dc061..6cd48a5 100644
--- a/odespy/demos/demo_Lsodi_1.py
+++ b/odespy/demos/demo_Lsodi_1.py
@@ -9,6 +9,7 @@
 This example is the typical usage of Lsodi with 
 user-supplied functions composed in Python.
 """
+from __future__ import print_function
 from odespy import *
 import scitools.std as st
 import numpy as np
@@ -52,7 +53,7 @@ def jac(u, t, s):
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'b-', title="Lsodi with Python functions",
         legend="with res, adda, ydoti & jac", hold="on")
-print 'Max error for test case 1 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(u[-1] - exact_final))
 
 # Test case 2: Lsodi, with res, ydoti & adda
 m = method(res=res, rtol=rtol, atol=atol, ydoti=ydoti,
@@ -61,6 +62,6 @@ def jac(u, t, s):
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'r*', title="Lsodi with Python functions",
         legend="with res, adda & ydoti", hold="on")
-print 'Max error for test case 1 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(u[-1] - exact_final))
 
 
diff --git a/odespy/demos/demo_Lsodi_1_fortran.py b/odespy/demos/demo_Lsodi_1_fortran.py
index c6d0203..ee0b592 100644
--- a/odespy/demos/demo_Lsodi_1_fortran.py
+++ b/odespy/demos/demo_Lsodi_1_fortran.py
@@ -8,6 +8,7 @@
 user-supplied functions composed in Fortran code.
 
 """
+from __future__ import print_function
 from odespy import *
 import scitools.std as st
 import numpy as np
@@ -83,7 +84,7 @@
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'r-', title="Lsodi with Fortran subroutines",
         legend="with res, adda, ydoti & jac", hold="on")
-print 'Max error for test case 1 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(u[-1] - exact_final))
 
 # Test case 2: Lsodi, with res, ydoti & adda
 m = method(res_f77=res_f77, rtol=rtol, atol=atol, ydoti=ydoti,
@@ -92,6 +93,6 @@
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'g*', title="Lsodi with Fortran subroutines",
         legend="with res, adda & ydoti", hold="on")
-print 'Max error for test case 2 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 2 is %g' % max(u[-1] - exact_final))
 
 os.remove('callback.so')
diff --git a/odespy/demos/demo_Lsodi_2.py b/odespy/demos/demo_Lsodi_2.py
index 2a92fb5..057816e 100644
--- a/odespy/demos/demo_Lsodi_2.py
+++ b/odespy/demos/demo_Lsodi_2.py
@@ -27,6 +27,7 @@
         0.5000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
 
 """
+from __future__ import print_function
 from odespy import *
 #import scitools.basics,easyviz as st
 import scitools.std as st
@@ -113,7 +114,7 @@ def jac_banded(u, t, s, ml, mu):   # Banded jacobian
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'r-', title="Lsodi with Python functions",
         legend="with res, adda_full & jac_full", hold="on")
-print 'Max error with test case 1 is %g' % max(u[-1] - exact_final)
+print('Max error with test case 1 is %g' % max(u[-1] - exact_final))
 
 
 # Test case 2: Lsodi, with res & adda_full
@@ -122,7 +123,7 @@ def jac_banded(u, t, s, ml, mu):   # Banded jacobian
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'b*', title="Lsodi with Python functions",
         legend="with res & adda_full", hold="on")
-print 'Max error with test case 2 is %g' % max(u[-1] - exact_final)
+print('Max error with test case 2 is %g' % max(u[-1] - exact_final))
 
 # Test case 3: Lsodi, with res, adda_banded, ml, mu, jac_banded
 m = method(res=res, rtol=rtol, atol=atol, 
@@ -133,7 +134,7 @@ def jac_banded(u, t, s, ml, mu):   # Banded jacobian
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'go', title="Lsodi with Python functions",
         legend="with res, adda_banded, jac_banded, ml, mu", hold="on")
-print 'Max error with test case 3 is %g' % max(u[-1] - exact_final)
+print('Max error with test case 3 is %g' % max(u[-1] - exact_final))
 
 # Test case 4: Lsodi, with res, adda_banded, ml, mu
 m = method(res=res, rtol=rtol, atol=atol, 
@@ -143,7 +144,7 @@ def jac_banded(u, t, s, ml, mu):   # Banded jacobian
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'y-', title="Lsodi with Python functions",
         legend="with res, adda_banded, ml, mu", hold="on")
-print 'Max error with test case 4 is %g' % max(u[-1] - exact_final)
+print('Max error with test case 4 is %g' % max(u[-1] - exact_final))
 
 
 
diff --git a/odespy/demos/demo_Lsodis_1.py b/odespy/demos/demo_Lsodis_1.py
index 9894c5c..6a93c67 100644
--- a/odespy/demos/demo_Lsodis_1.py
+++ b/odespy/demos/demo_Lsodis_1.py
@@ -26,6 +26,7 @@
      = r4d*(y(nm1)**2-y(1)**2)+eodsq*(y(1)-2*y(NEQ)+y(nm1))
      where r4d = 1/(4*delx), eodsq = eta/delx**2 and nm1 = NEQ-1.
 """
+from __future__ import print_function
 from odespy import *
 import numpy as np
 
@@ -73,11 +74,11 @@ def jac(y, t, s, j, ia, ja):
 m = method(res=res, adda_lsodis=adda, atol=atol, rtol=rtol, jac_lsodis=jac)
 m.set_initial_condition(u0)
 y, t = m.solve(time_points)
-print 'Max error for test case 1 is %g' % max(y[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(y[-1] - exact_final))
 
 # Test case 2:  With res & adda
 m = method(res=res, adda_lsodis=adda, atol=atol, rtol=rtol,lrw=4000,liw=100)
 m.set_initial_condition(u0)
 y, t = m.solve(time_points)
-print 'Max error for test case 2 is %g' % max(y[-1] - exact_final)
+print('Max error for test case 2 is %g' % max(y[-1] - exact_final))
 
diff --git a/odespy/demos/demo_Lsodis_2.py b/odespy/demos/demo_Lsodis_2.py
index aa5adfb..7c3a649 100644
--- a/odespy/demos/demo_Lsodis_2.py
+++ b/odespy/demos/demo_Lsodis_2.py
@@ -19,6 +19,7 @@
     An ODE system is generated by a simplified Galerkin treatment of the spatial
     variable x.
 """
+from __future__ import print_function
 from odespy import *
 import numpy as np
 
@@ -69,25 +70,25 @@ def jac(y, t, s, j, ia, ja):
            ia=ia, ja=ja, ic=ic, jc=jc)
 m.set_initial_condition(u0)
 y, t = m.solve(time_points)
-print 'Max error for test case 1 is %g' % max(y[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(y[-1] - exact_final))
 
 # Test case 2: with res, adda, ia, ja & jac
 m = method(res=res, adda_lsodis=adda, atol=atol, rtol=rtol, jac_lsodis=jac,
            ia=ia, ja=ja)
 m.set_initial_condition(u0)
 y, t = m.solve(time_points)
-print 'Max error for test case 2 is %g' % max(y[-1] - exact_final)
+print('Max error for test case 2 is %g' % max(y[-1] - exact_final))
 
 # Test case 3: with res, adda & jac
 m = method(res=res, adda_lsodis=adda, atol=atol, rtol=rtol, jac_lsodis=jac)
 m.set_initial_condition(u0)
 y, t = m.solve(time_points)
-print 'Max error for test case 3 is %g' % max(y[-1] - exact_final)
+print('Max error for test case 3 is %g' % max(y[-1] - exact_final))
 
 # Test case 4:  With res & adda
 m = method(res=res, adda_lsodis=adda, atol=atol, rtol=rtol,lrw=4000,liw=100)
 m.set_initial_condition(u0)
 y, t = m.solve(time_points)
-print 'Max error for test case 4 is %g' % max(y[-1] - exact_final)
+print('Max error for test case 4 is %g' % max(y[-1] - exact_final))
 
 
diff --git a/odespy/demos/demo_Lsoibt_2.py b/odespy/demos/demo_Lsoibt_2.py
index e593dd5..0dfbb3b 100644
--- a/odespy/demos/demo_Lsoibt_2.py
+++ b/odespy/demos/demo_Lsoibt_2.py
@@ -32,6 +32,7 @@
 Terms involving boundary values (subscripts 0 or 100) are dropped
 from the equations for k = 1 and k = 99 above.
 """
+from __future__ import print_function
 from odespy import *
 import scitools.std as st
 import numpy as np
@@ -209,7 +210,7 @@ def jac(y, t, s):
 u_final = u[-1].reshape(99,3)
 u1, u2, u3 = u_final[:, 0], u_final[:, 1], u_final[:, 2]
 max_error = max(max(u1 - u1_exact), max(u2 - u2_exact), max(u3 - u3_exact))
-print 'Max error with Test case 1 is %g' % max_error
+print('Max error with Test case 1 is %g' % max_error)
 
 # Test case 2: Lsoibt, with res, adda, mb, nb
 m = method(rtol=rtol, atol=atol, res=res, adda_lsoibt=adda, 
@@ -222,4 +223,4 @@ def jac(y, t, s):
 u_final = u[-1].reshape(99,3)
 u1, u2, u3 = u_final[:, 0], u_final[:, 1], u_final[:, 2]
 max_error = max(max(u1 - u1_exact), max(u2 - u2_exact), max(u3 - u3_exact))
-print 'Max error with Test case 2 is %g' % max_error
+print('Max error with Test case 2 is %g' % max_error)
diff --git a/odespy/demos/demo_MyRungeKutta.py b/odespy/demos/demo_MyRungeKutta.py
index 6353bc5..dd35918 100644
--- a/odespy/demos/demo_MyRungeKutta.py
+++ b/odespy/demos/demo_MyRungeKutta.py
@@ -4,6 +4,7 @@
 This example intends to show users how to apply MyRungeKutta to define 
 own RungKutta solvers.
 """
+from __future__ import print_function
 from odespy import *
 #import scitools.basics,easyviz as st
 import scitools.std as st
@@ -45,7 +46,7 @@ def f(u,t):
     method.set_initial_condition(u0)
     u,t = method.solve(time_points)
     error = abs((u[-1] - np.exp(-1.))/np.exp(-1.))
-    print 'Error is %g with solver %s' % (error, m)
+    print('Error is %g with solver %s' % (error, m))
 
 
 
diff --git a/odespy/demos/demo_RKC_2.py b/odespy/demos/demo_RKC_2.py
index 0cfcac3..7cf6914 100644
--- a/odespy/demos/demo_RKC_2.py
+++ b/odespy/demos/demo_RKC_2.py
@@ -14,6 +14,7 @@
   in Fortran language instead of Python. See demo_RKC_2_fortran.py.
   
 """
+from __future__ import print_function
 
 from odespy import *
 import scitools.std as st
@@ -74,7 +75,7 @@ def spcrad(u,t):
 
 st.figure()
 method = RKC
-print 'This test will take several minutes, please wait...'
+print('This test will take several minutes, please wait...')
 
 # Test case 1: RKC with f, jac_constant, spcrad
 m = method(f, rtol=rtol, atol=atol, spcrad=spcrad, jac_constant=jac_constant)
diff --git a/odespy/demos/demo_Radau5_1.py b/odespy/demos/demo_Radau5_1.py
index 3b6b943..dd8cfa1 100644
--- a/odespy/demos/demo_Radau5_1.py
+++ b/odespy/demos/demo_Radau5_1.py
@@ -1,3 +1,4 @@
+from __future__ import print_function
 from odespy import *
 import scitools.std as st
 
@@ -29,7 +30,7 @@ def jac(u,t):
 u0[0], u0[-1] = 1., .0057
 time_points = np.linspace(0., 10., 20)
 
-print "HIRES, chemical reaction, mildly stiff"
+print("HIRES, chemical reaction, mildly stiff")
 st.figure()
 
 # Loop for all possible solvers
@@ -38,4 +39,4 @@ def jac(u,t):
     method.set_initial_condition(u0)
     u,t = method.solve(time_points)
     st.plot(t, u[:,0], hold="on", legend=solver, axis=[0.,10.,0.,1.])
-    print 'Succeed when solver is %s' % solver
+    print('Succeed when solver is %s' % solver)
diff --git a/odespy/demos/demo_Radau5_4.py b/odespy/demos/demo_Radau5_4.py
index d068337..2b3dc4f 100644
--- a/odespy/demos/demo_Radau5_4.py
+++ b/odespy/demos/demo_Radau5_4.py
@@ -9,6 +9,7 @@
 This example is the typical usage of DAE with
 user-supplied functions composed in Python.
 """
+from __future__ import print_function
 from odespy import *
 import scitools.std as st
 import numpy as np
@@ -51,7 +52,7 @@ def jac(u, t):
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'b-', title="Radau5 with Python functions",
         legend="with f, mas & jac", hold="on")
-print 'Max error for test case 1 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(u[-1] - exact_final))
 
 # Test case 2: Radau5, with f, mas
 m = method(f=f, mas=mas, rtol=rtol, atol=atol)
@@ -59,6 +60,6 @@ def jac(u, t):
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'r*', title="Radau5 with Python functions",
         legend="with f, mas", hold="on")
-print 'Max error for test case 1 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(u[-1] - exact_final))
 
 
diff --git a/odespy/demos/demo_Radau5_4_fortran.py b/odespy/demos/demo_Radau5_4_fortran.py
index 2c40db3..8b36b66 100644
--- a/odespy/demos/demo_Radau5_4_fortran.py
+++ b/odespy/demos/demo_Radau5_4_fortran.py
@@ -8,6 +8,7 @@
 user-supplied functions composed in Fortran code.
 
 """
+from __future__ import print_function
 from odespy import *
 import scitools.std as st
 import numpy as np
@@ -91,7 +92,7 @@
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'r-', title="Radau5 with Fortran subroutines",
         legend="with f, mas & jac", hold="on")
-print 'Max error for test case 1 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(u[-1] - exact_final))
 
 # Test case 2: Radau5, with f & mas
 m = method(None, f_f77=f_str, rtol=rtol, atol=atol, mas_f77=mas_str)
@@ -99,6 +100,6 @@
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'g*', title="Radau5 with Fortran subroutines",
         legend="with f & mas", hold="on")
-print 'Max error for test case 2 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 2 is %g' % max(u[-1] - exact_final))
 
 os.remove('tmp_callback.so')
diff --git a/odespy/demos/demo_Radau5_5.py b/odespy/demos/demo_Radau5_5.py
index 2476975..c29fa0f 100644
--- a/odespy/demos/demo_Radau5_5.py
+++ b/odespy/demos/demo_Radau5_5.py
@@ -26,6 +26,7 @@
         0.0000E+00  0.0000E+00  0.5000E+00  0.1000E+01  0.1000E+01  0.1000E+01
         0.5000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00
 """
+from __future__ import print_function
 
 from odespy import *
 import scitools.std as st
@@ -96,7 +97,7 @@ def jac_banded(u, t, ml, mu):   # Banded jacobian
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'b-', title="Radau5 with Python functions",
         legend="with f, mas & jac_full", hold="on")
-print 'Max error for test case 1 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(u[-1] - exact_final))
 
 # Test case 2: Radau5, with f, mas
 m = method(f=f, mas=mas, rtol=rtol, atol=atol)
@@ -104,7 +105,7 @@ def jac_banded(u, t, ml, mu):   # Banded jacobian
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'r*', title="Radau5 with Python functions",
         legend="with f, mas", hold="on")
-print 'Max error for test case 2 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 2 is %g' % max(u[-1] - exact_final))
 
 # Test case 3: Radau5, with f, mas_banded, ml, mu & jac_banded
 m = method(f=f, mas=mas_banded, rtol=rtol, atol=atol,
@@ -114,6 +115,6 @@ def jac_banded(u, t, ml, mu):   # Banded jacobian
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'b-', title="Radau5 with Python functions",
         legend="with f, mas & jac_banded", hold="on")
-print 'Max error for test case 3 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 3 is %g' % max(u[-1] - exact_final))
 
 
diff --git a/odespy/demos/demo_Radau5_5_fortran.py b/odespy/demos/demo_Radau5_5_fortran.py
index 7aea9e9..6dc3ee0 100644
--- a/odespy/demos/demo_Radau5_5_fortran.py
+++ b/odespy/demos/demo_Radau5_5_fortran.py
@@ -8,6 +8,7 @@
 user-supplied functions composed in Fortran code.
 
 """
+from __future__ import print_function
 from odespy import *
 import scitools.std as st
 import numpy as np
@@ -141,7 +142,7 @@
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'b-', title="Radau5 with Fortran subroutines",
         legend="with f, mas & jac_full", hold="on")
-print 'Max error for test case 1 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 1 is %g' % max(u[-1] - exact_final))
 
 # Test case 2: Radau5, with f, mas
 m = method(None, f_f77=f_f77, mas_f77=mas_f77, rtol=rtol, atol=atol)
@@ -149,7 +150,7 @@
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'r*', title="Radau5 with Fortran subroutines",
         legend="with f, mas", hold="on")
-print 'Max error for test case 2 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 2 is %g' % max(u[-1] - exact_final))
 
 # compile these Fortran subroutines
 os.remove('tmp_callback.so')
@@ -170,7 +171,7 @@
 u,t = m.solve(time_points)
 st.plot(t, u[:,0], 'b-', title="Radau5 with Fortran subroutines",
         legend="with f, mas_banded & jac_banded", hold="on")
-print 'Max error for test case 3 is %g' % max(u[-1] - exact_final)
+print('Max error for test case 3 is %g' % max(u[-1] - exact_final))
 
 os.remove('tmp_callback2.so')
 
diff --git a/odespy/demos/demo_basics_complex_Pi.py b/odespy/demos/demo_basics_complex_Pi.py
index 0b91200..c7cdf91 100644
--- a/odespy/demos/demo_basics_complex_Pi.py
+++ b/odespy/demos/demo_basics_complex_Pi.py
@@ -1,3 +1,4 @@
+from __future__ import print_function
 # Author: Liwei Wang
 # Scalar ODE with complex value
 # u' = 1/(t - 10 + 1j)
@@ -18,8 +19,8 @@
 f = lambda u, t: 1./(t - 10. + 1j)
 u0, time_points = 0, np.linspace(0., 20., 200)
 
-print """ Scalar ODE with complex value
- u' =  1./(t - 10 + 1j) """
+print(""" Scalar ODE with complex value
+ u' =  1./(t - 10 + 1j) """)
 
 succeed_no = 0
 for solver in complex_solvers:
@@ -28,6 +29,6 @@
         method.set_initial_condition(u0)
         u,t = method.solve(time_points)
         success = np.allclose(u[-1].imag, -2*np.arctan(10))
-        print solver, 'Succeed' if success else 'Fail'
+        print(solver, 'Succeed' if success else 'Fail')
     except:
-        print 'Failed when solver is %s' % solver
+        print('Failed when solver is %s' % solver)
diff --git a/odespy/demos/demo_basics_exponential.py b/odespy/demos/demo_basics_exponential.py
index 7354789..1aa0b7e 100644
--- a/odespy/demos/demo_basics_exponential.py
+++ b/odespy/demos/demo_basics_exponential.py
@@ -1,3 +1,4 @@
+from __future__ import print_function
 # Author: Liwei Wang
 
 # Scalar ODE:  Exponential
@@ -17,7 +18,7 @@
 f=lambda u,t:-u
 u0, time_points = 1., np.linspace(0., 10., 100)
 
-print """Scalar ODE: Exponential u = exp(-t), u' = -u"""
+print("""Scalar ODE: Exponential u = exp(-t), u' = -u""")
 
 # Loop for all possible solvers
 for solver in solvers:
@@ -30,6 +31,6 @@
             st.figure()           # Initialize new figure.
         st.plot(t, u, hold="on", legend=solver, axis=[0., 10., 0., 1.5])
         solver_no += 1
-        print 'Succeed when solver is %s' % solver
+        print('Succeed when solver is %s' % solver)
     except:
-        print 'Failed when solver is %s' % solver
+        print('Failed when solver is %s' % solver)
diff --git a/odespy/demos/demo_basics_sine.py b/odespy/demos/demo_basics_sine.py
index 25b1c7e..587befa 100644
--- a/odespy/demos/demo_basics_sine.py
+++ b/odespy/demos/demo_basics_sine.py
@@ -1,3 +1,4 @@
+from __future__ import print_function
 # Nonstiff ODE:  Sine
 # u'' = - u,   u = sin(t)
 
@@ -17,7 +18,7 @@
 
 u0, time_points = [0., 1.], np.linspace(0., 10., 100)
 
-print """Nonstiff ODE: Sine u = sin(t), u'' = -u"""
+print("""Nonstiff ODE: Sine u = sin(t), u'' = -u""")
 
 # Loop for all possible solvers
 for solver in solvers:
@@ -30,7 +31,7 @@
             solver_no = 1
         st.plot(t, u[:,0], hold="on", legend=solver, axis=[0., 10., -1.5, 1.5])
         solver_no += 1
-        print 'Succeed when solver is %s' % solver
+        print('Succeed when solver is %s' % solver)
     except:
-        print 'Failed when solver is %s' % solver
+        print('Failed when solver is %s' % solver)
 
diff --git a/odespy/demos/demo_basics_van_der_pol.py b/odespy/demos/demo_basics_van_der_pol.py
index c1ac32a..ad10cf0 100644
--- a/odespy/demos/demo_basics_van_der_pol.py
+++ b/odespy/demos/demo_basics_van_der_pol.py
@@ -1,3 +1,4 @@
+from __future__ import print_function
 # Stiff ODE:  Van der Pol oscillator
 # u'' = 3*(1 - u**2)*u' - u
 from odespy import *
@@ -16,8 +17,8 @@
 
 u0, time_points = [2., 0.], np.linspace(0., 10., 100)
 
-print """Van der Pol oscillator problem:
-     u'' = 3*(1 - u**2)*u' - u"""
+print("""Van der Pol oscillator problem:
+     u'' = 3*(1 - u**2)*u' - u""")
 
 # Loop for all possible solvers
 for solver in solvers:
@@ -30,6 +31,6 @@
             solver_no = 1
         st.plot(t, u[:,0], hold="on", legend=solver, axis=[0.,10.,-4.,4.])
         solver_no += 1
-        print 'Succeed when solver is %s' % solver
+        print('Succeed when solver is %s' % solver)
     except:
-        print 'Failed when solver is %s' % solver
+        print('Failed when solver is %s' % solver)
diff --git a/odespy/demos/demo_odelab.py b/odespy/demos/demo_odelab.py
index 6efc74c..44056c1 100644
--- a/odespy/demos/demo_odelab.py
+++ b/odespy/demos/demo_odelab.py
@@ -1,3 +1,4 @@
+from __future__ import print_function
 import odespy, numpy, matplotlib.pyplot as plt
 
 # Run through all scheme classes in with a simple ODE u'=-u, u(0)=1
@@ -11,18 +12,18 @@
 
 failures = []
 for scheme in schemes:
-    print scheme,
+    print(scheme, end=' ')
     solver = odespy.odelab(problem.f, odelab_solver=scheme,
                            atol=1E-2, rtol=1E-2)
     solver.set_initial_condition(1.0)
     try:
         u, t = solver.solve(numpy.linspace(0, 3, 221))
-        print 'summary:', len(u), u[-1], numpy.exp(-t[-1]) - u[-1]
+        print('summary:', len(u), u[-1], numpy.exp(-t[-1]) - u[-1])
     except:
         failures.append(scheme)
 
-print '\nFailed classes:\n', ' '.join(failures)
-print '\nAll solver classes:\n', ' '.join(schemes)
-print '\nWorking classes:\n', ' '.join([scheme for scheme in schemes if scheme not in failures])
-print '\nDocumented working classes:\n', ' '.join(odespy.odelab.solvers)
+print('\nFailed classes:\n', ' '.join(failures))
+print('\nAll solver classes:\n', ' '.join(schemes))
+print('\nWorking classes:\n', ' '.join([scheme for scheme in schemes if scheme not in failures]))
+print('\nDocumented working classes:\n', ' '.join(odespy.odelab.solvers))
 
diff --git a/odespy/odepack.py b/odespy/odepack.py
index 96523c6..95a4d19 100644
--- a/odespy/odepack.py
+++ b/odespy/odepack.py
@@ -1,6 +1,8 @@
+from __future__ import print_function
+from __future__ import absolute_import
 __author__ = ['Liwei Wang, Univ. of Oslo']
 
-from solvers import Solver
+from .solvers import Solver
 import numpy as np
 import sys, inspect
 
@@ -463,7 +465,7 @@
 
     )
 
-import solvers
+from . import solvers
 solvers._parameters.update(_parameters_Odepack)
 
 
@@ -735,9 +737,9 @@ def check_liwlrw(self):
             else:
                 value = getattr(self, name)
                 if value < min_value:
-                    print '''
+                    print('''
           Insufficient input! "%s"=%d are reset to be the minimum size = %d '''\
-                        % (name, value, min_value)
+                        % (name, value, min_value))
                     setattr(self, name, min_value)
 
     def check_tol(self):
@@ -759,8 +761,8 @@ def check_pars(self):
         for pars in (('ia', 'ja'), ('ic', 'jc'), ('ml', 'mu'), ('mb', 'nb')):
             arg_a, arg_b = pars
             if int(hasattr(self, arg_a) + hasattr(self, arg_b)) == 1:
-                raise ValueError,'\
-        Error! %s and %s have to be input simutaneously!' % (arg_a, arg_b)
+                raise ValueError('\
+        Error! %s and %s have to be input simutaneously!' % (arg_a, arg_b))
 
     def check_iaja(self):
         '''
@@ -797,7 +799,7 @@ def check_iaja(self):
                               array_a[-1] != len(array_b)+1)
                 for error_index in range(3):
                     if iaja_check[error_index]:
-                        raise ValueError, err_messages[error_index]
+                        raise ValueError(err_messages[error_index])
 
     def validate_data(self):
         '''
@@ -945,17 +947,17 @@ def new_stepnr(self):
         '''
         nsteps = getattr(self, 'nsteps', 500)
         if nsteps == 2000:    # The maximum step-amount has been reached
-            raise ValueError, '''
+            raise ValueError('''
         Failed iteration although step number has been set to 2000.
-        Please check your input.'''
+        Please check your input.''')
         mx_new = min(nsteps+200, 2000)
         # maximum limitation is set to 2000
         self.nsteps = mx_new  # valid for the following steps
-        print '''\
+        print('''\
         Excessive amount of work detected!
         Input step amount "nsteps"=%d is not enough for iteration!
         nsteps has been reset to %d to avoid this error!'''\
-        % (nsteps, mx_new)
+        % (nsteps, mx_new))
         return mx_new
 
     def tol_multiply(self, tolsf):
@@ -966,16 +968,16 @@ def tol_multiply(self, tolsf):
         Then we could try to adjust tolerance settings with suggested factor
         to avoid this error.
         '''
-        print 'Tolerance is scaled by suggested factor %.2g''' % tolsf
+        print('Tolerance is scaled by suggested factor %.2g''' % tolsf)
         self.rtol *= tolsf
         self.atol *= tolsf
 
     def expand_rwork(self, new_lrw, expand=False):
         '''
         Length of real work array is smaller than actually required length.
-	Then we could expand work array to avoid this error.
+        Then we could expand work array to avoid this error.
         '''
-        print 'The length of real work array has been reset to %d' % new_lrw
+        print('The length of real work array has been reset to %d' % new_lrw)
         if expand:          # Expand real arrays for linearly implicit solvers
             self.rwork = list(self.rwork) + [0.]*(new_lrw-self.lrw)
         self.lrw = new_lrw
@@ -984,9 +986,9 @@ def expand_iwork(self, new_liw, expand=False):
         '''
         Extension module return an actually required length for
         integer work array when it is too short.
-	Then we could expand work array to required length to avoid this error.
+        Then we could expand work array to required length to avoid this error.
         '''
-        print 'The length of integer work array has been reset to %d' % new_liw
+        print('The length of integer work array has been reset to %d' % new_liw)
         if expand:          # Expand integer arrays for linearly implicit solvers
             self.iwork = list(self.iwork) + [0.]*(new_liw - self.liw)
         self.liw = new_liw
@@ -1023,11 +1025,11 @@ def print_roots(self, jroot, t_current, u_current):
             value = g(t_current, u_current)
         for i in range(ng):
             if jroot[i]:   # found root for i-th constraint equation
-                print '''
+                print('''
         Root found at t = %g for %dth constraint function in g''' \
-                    % (t_current, i+1)
+                    % (t_current, i+1))
                 value_ith = value[i] if ng > 1 else value
-                print 'Error in location of root is %g' % value_ith
+                print('Error in location of root is %g' % value_ith)
 
 
     def solve(self, time_points, terminate=None):
@@ -1069,13 +1071,11 @@ def solve(self, time_points, terminate=None):
                                  self.jac_column_f77)
             # call extension module
             nstop, u, istate, rinfo, iinfo = \
-                apply(self._odepack.solve,
-                      (terminate_int, itermin, nstop, f, u,
+                self._odepack.solve(*(terminate_int, itermin, nstop, f, u,
                        self.t, self.itol, self.rtol, self.atol,
                        istate, self.iopt, self.rwork_in, self.lrw,
                        self.iwork_in, self.liw, jac, jac_column,
-                       self.mf, g, self.ng, solver_name),
-                      self._extra_args_fortran)
+                       self.mf, g, self.ng, solver_name), **self._extra_args_fortran)
             tried += 1
             if nstop == step_no:    # successful
                 self.finished = True
@@ -1087,11 +1087,11 @@ def solve(self, time_points, terminate=None):
                 self.finished = True
                 tried = 0
             elif istate == 0:
-                print "Iteration stops at step Nr.%d," % nstop
-                print " when function TERMINATE return with True."
+                print("Iteration stops at step Nr.%d," % nstop)
+                print(" when function TERMINATE return with True.")
                 self.finished = True
             elif istate < 0:                   # Error occurs!
-                print self._error_messages[istate] + str(rinfo[1])
+                print(self._error_messages[istate] + str(rinfo[1]))
                 if istate == -1:    # Increase maximum step-number.
                     self.iwork_in[5] = self.new_stepnr()
                     self.iopt = 1
@@ -1168,19 +1168,16 @@ def advance(self):
 
         tried = 0
         while tried < 5:       # prevent endless loop
-            u_new, t, istate, iwork = apply(\
-                eval('self._odepack.d%s' % solver_name),\
-                (res, adda, jac, neq, u[n].copy(), self.ydoti, t, t_next,
+            u_new, t, istate, iwork = eval('self._odepack.d%s' % solver_name)(*(res, adda, jac, neq, u[n].copy(), self.ydoti, t, t_next,
                  self.itol, self.rtol, self.atol, itask, istate, self.iopt,
-                 self.rwork, self.lrw, self.iwork, self.liw, self.mf),
-                self._extra_args_fortran)
+                 self.rwork, self.lrw, self.iwork, self.liw, self.mf), **self._extra_args_fortran)
             tried += 1
-	    # "istate" indicates the returned status
+            # "istate" indicates the returned status
             if istate >= 1:
                 # successful return status
                 break
             else:       # Error occurs!
-                print self._error_messages[istate] + str(self.rwork[12])
+                print(self._error_messages[istate] + str(self.rwork[12]))
                 if istate == -1:    # Increase maximum step-number.
                     self.iwork[5], self.iopt = self.new_stepnr(), 1
                 elif istate == -2:  # Multiply tolerance with suggested factor
@@ -1200,7 +1197,7 @@ def advance(self):
                 else:   # Unavoidable interrupts
                     sys.exit(1)  #  Interrupt
                 istate = 1
-	return u_new
+        return u_new
 
 ### end of class Odepack ###
 
@@ -1229,7 +1226,7 @@ class Lsode(Odepack):
     _extra_args_fortran = {}
 
     def adjust_parameters(self):
-	"""Properties for new parameters in this solver."""
+        """Properties for new parameters in this solver."""
        # If jac_banded is input in form of jac(u,t,ml,mu),
         # wrap jac_banded to jac_banded_f77 for Fortran code
         self._parameters['jac_banded']['paralist_old'] = 'u,t,ml,mu'
@@ -1267,7 +1264,7 @@ def adjust_parameters(self):
         Odepack.adjust_parameters(self)
 
     def set_extra_args(self):
-	# ml & mu are required to be extra parameters for banded Jacobian.
+        # ml & mu are required to be extra parameters for banded Jacobian.
         if hasattr(self,'ml') and hasattr(self,'mu'):
             if self.iter_method == 4 and \
                     (not hasattr(self.jac_f77, '_cpointer')):
@@ -1340,7 +1337,7 @@ class Lsoda(Odepack):
     was successful as far as '''
 
     def adjust_parameters(self):
-	"""Properties for new parameters in this solver."""
+        """Properties for new parameters in this solver."""
         # If jac_banded is input in form of jac(u,t,ml,mu),
         # wrap jac_banded to jac_banded_f77 for Fortran code
         self._parameters['jac_banded']['paralist_old'] = 'u,t,ml,mu'
@@ -1369,7 +1366,7 @@ def adjust_parameters(self):
         Odepack.adjust_parameters(self)
 
     def set_extra_args(self):
-	# ml & mu are required to be extra parameters for banded Jacobian.
+        # ml & mu are required to be extra parameters for banded Jacobian.
         if hasattr(self,'ml') and hasattr(self,'mu'):
             if self.iter_method == 4 and \
                     (not hasattr(self.jac_f77, '_cpointer')):
@@ -1453,7 +1450,7 @@ class Lsodar(Odepack):
     was successful as far as '''
 
     def adjust_parameters(self):
-	"""Properties for new parameters in this solver."""
+        """Properties for new parameters in this solver."""
         # If jac_banded is input in form of jac(u,t,ml,mu),
         # wrap jac_banded to jac_banded_f77 for Fortran code
         self._parameters['jac_banded']['paralist_old'] = 'u,t,ml,mu'
@@ -1494,7 +1491,7 @@ def adjust_parameters(self):
 
 
     def set_extra_args(self):
-	# ml & mu are required to be extra parameters for banded Jacobian.
+        # ml & mu are required to be extra parameters for banded Jacobian.
         if hasattr(self,'ml') and hasattr(self,'mu'):
             if self.iter_method == 4 and \
                     (not hasattr(self.jac_f77, '_cpointer')):
@@ -1541,9 +1538,9 @@ def solve(self, time_points, terminate=None):
             if not hasattr(self.g_f77,'_cpointer'):
                 self.ng = np.asarray(self.g(time_points[0],self.U0)).size
             elif not hasattr(self,'ng'):
-                raise ValueError, '''
+                raise ValueError('''
         Unsufficient input! ng must be specified if g is input as a
-        Fortran subroutine. '''
+        Fortran subroutine. ''')
         return Odepack.solve(self,time_points, terminate=terminate)
 
 ### End of Lsodar ###
@@ -1571,7 +1568,7 @@ class Lsodes(Odepack):
 
 
     def adjust_parameters(self):
-	"""Properties for new parameters in this solver."""
+        """Properties for new parameters in this solver."""
         # If jac_column is input in form of jac(u,t,j),
         # wrap jac_column to jac_column_f77(t,u,j-1) for Fortran code.
         self._parameters['jac_column']['paralist_old'] = 'u,t,j-1,ia,ja'
@@ -1777,7 +1774,7 @@ def set_iter_method(self):
             elif with_full_adda:
                 self.iter_method = 2
             else:
-                raise ValueError, 'adda must be supplied in Lsodi.'
+                raise ValueError('adda must be supplied in Lsodi.')
 
     def set_jac(self):
         if self.iter_method == 4:
@@ -2072,9 +2069,9 @@ def validate_data(self):
         if not Odepack.validate_data(self):
             return False
         if self.mb*self.nb != self.neq:
-                raise ValueError,'''
+                raise ValueError('''
     The requirement for block size (mb,nb) are: mb>=1, nb>=4, mb*nb=neq=%d.
-    Your block size are (%d,%d). ''' % (self.neq, mb, nb)
+    Your block size are (%d,%d). ''' % (self.neq, mb, nb))
         return True
 
     def solve(self, time_points, terminate=None):
diff --git a/odespy/problems.py b/odespy/problems.py
index 023e218..740ac53 100644
--- a/odespy/problems.py
+++ b/odespy/problems.py
@@ -1,5 +1,7 @@
+from __future__ import print_function
+from __future__ import absolute_import
 import numpy as np
-from solvers import compile_f77
+from .solvers import compile_f77
 
 class Problem:
     '''
@@ -282,7 +284,7 @@ def u_exact(self, t):
 def default_oscillator(P, resolution_per_period=20):
     n = 4.5
     r = resolution_per_period # short form
-    print 'default', P, n*P
+    print('default', P, n*P)
     tp = np.linspace(0, n*P, r*n+1)
     # atol=rtol since u approx 1, atol level set at the
     # error RK4 produces with 20 steps per period, times 0.05
@@ -519,7 +521,7 @@ def tester(problems, methods, time_points=None, compare_tol=1E-4,
     error_msg = {}
     for problem in problems:
         pname = problem.__class__.__name__
-        print 'problem ', pname
+        print('problem ', pname)
         methods4problem = [method for method in methods
                            if method not in problem.not_suitable_solvers]
         defaults = problem.default_parameters()
@@ -541,7 +543,7 @@ def tester(problems, methods, time_points=None, compare_tol=1E-4,
                 this_solver_prm[name] = defaults[name]
 
         for method in methods4problem:
-            print '  testing', method,
+            print('  testing', method, end=' ')
             solver = eval('odespy.'+method)(problem.f)
 
             # Important to set parameters before setting initial cond.
@@ -554,12 +556,12 @@ def tester(problems, methods, time_points=None, compare_tol=1E-4,
             error = problem.verify(u, t)
             if error is not None:
                 results[pname][method] = (u, error)
-                print error,
+                print(error, end=' ')
                 if error > compare_tol:
-                    print 'WARNING: tolerance %.0E exceeded' % compare_tol,
+                    print('WARNING: tolerance %.0E exceeded' % compare_tol, end=' ')
             else:
                 results[pname][method] = (u,)
-            print
+            print()
     return results
 
 """
diff --git a/odespy/problems_pde.py b/odespy/problems_pde.py
index ce9b266..c0c245d 100644
--- a/odespy/problems_pde.py
+++ b/odespy/problems_pde.py
@@ -1,4 +1,5 @@
-from problems import Problem, np
+from __future__ import absolute_import
+from .problems import Problem, np
 
 class Diffusion1D(Problem):
     """
diff --git a/odespy/radau5.py b/odespy/radau5.py
index f185888..530b986 100644
--- a/odespy/radau5.py
+++ b/odespy/radau5.py
@@ -1,3 +1,5 @@
+from __future__ import print_function
+from __future__ import absolute_import
 __author__ = ['Liwei Wang, Univ. of Oslo']
 
 from odespy import Solver
@@ -80,7 +82,7 @@
         type=callable),
 )
 
-import solvers
+from . import solvers
 solvers._parameters.update(_parameters_Radau5)
 
 
@@ -278,9 +280,9 @@ def advance(self):
                              *args, **self._extra_args_fortran)
 
         if idid < 0:          # Error occurred
-            print self._error_messages[idid] + str(t_new)
+            print(self._error_messages[idid] + str(t_new))
             sys.exit(1)   # Interrupt
-	return u_new
+        return u_new
 
 ### end of class Radau5 ###
 
diff --git a/odespy/rkc.py b/odespy/rkc.py
index 92fb753..23740e9 100644
--- a/odespy/rkc.py
+++ b/odespy/rkc.py
@@ -1,6 +1,7 @@
 """Module for wrapping rkc.f."""
+from __future__ import absolute_import
 
-from solvers import Solver, Adaptive
+from .solvers import Solver, Adaptive
 import numpy as np
 
 # f_f77 and other items are defined in odepack.py and will
@@ -31,7 +32,7 @@
         type=callable),
     )
 
-import solvers
+from . import solvers
 solvers._parameters.update(_parameters_RKC)
 
 class RKC(Adaptive):
@@ -105,9 +106,9 @@ def check_atol(self):
         atol = self.atol
         if not isinstance(atol,float):
             if len(atol) not in (1, self.neq):
-                raise ValueError,  '''
+                raise ValueError('''
 ATOL =%s should be either a scalar or a vector of length NEQ=%d.
-           ''' % (str(atol), self.neq)
+           ''' % (str(atol), self.neq))
 
     def validate_data(self):
         self.check_atol()
diff --git a/odespy/rkf45.py b/odespy/rkf45.py
index 486ca30..c2d70a4 100644
--- a/odespy/rkf45.py
+++ b/odespy/rkf45.py
@@ -1,6 +1,7 @@
 """Module for wrapping rkf45.py."""
+from __future__ import absolute_import
 
-from solvers import Solver, Adaptive
+from .solvers import Solver, Adaptive
 import numpy as np
 
 
diff --git a/odespy/runtests.py b/odespy/runtests.py
index 7eb921c..e1a82c7 100644
--- a/odespy/runtests.py
+++ b/odespy/runtests.py
@@ -1,10 +1,12 @@
+from __future__ import print_function
+from __future__ import absolute_import
 from odespy import *
-from problems import *
+from .problems import *
 from scitools.std import plot, hold
 try:
     import joblib  # need joblib for storing database nested dict
 except ImportError:
-    print 'This script requires joblib to be installed.'
+    print('This script requires joblib to be installed.')
     sys.exit(1)
 
 newresults = {}
@@ -32,20 +34,20 @@ def test_disk():
     tp = np.linspace(0, 10, 1000000)
     # tp = problem.default_parameters()['time_points']
     u, t = solver.solve(tp)
-    print 'made u', u.size
+    print('made u', u.size)
     u0 = u[:,0]
-    print 'made u0'
+    print('made u0')
     if tp.size > 500000:
         u = u[-100:]
         t = t[-100:]
         error = problem.verify(u, t)
-        print error
+        print(error)
     else:
         error = None
-    print 'u memmap:', u[:10]
+    print('u memmap:', u[:10])
     u[:] = 5.5
     u.flush()
-    print u[:10]
+    print(u[:10])
     return {'Linear2': {'t': t, 'ForwardEuler': (u, error)}}
 
 def test_Exponential():
@@ -63,7 +65,7 @@ def test_VanDerPol(mu=0):
     methods = [Vode, RK4, RungeKutta4, ForwardEuler, BackwardEuler]
     for method in methods:
         name = method.__name__
-        print name
+        print(name)
         solver = method(problem.f, jac=problem.jac,
                         atol=d['atol'], rtol=d['rtol'])
         solver.set_initial_condition(problem.U0)
@@ -72,7 +74,7 @@ def test_VanDerPol(mu=0):
              legend_fancybox=True, legend_loc='upper left')
         hold('on')
         e = problem.verify(u, t)
-        if e is not None: print e
+        if e is not None: print(e)
 
 
 if __name__ == '__main__':
diff --git a/odespy/solvers.py b/odespy/solvers.py
index af595d4..dd34108 100644
--- a/odespy/solvers.py
+++ b/odespy/solvers.py
@@ -112,6 +112,7 @@ class ``BackwardEuler``, for instance.
 for doing one (adaptive) time step, called in the ``advance`` method.
 
 '''
+from __future__ import print_function
 
 import pprint, sys, os, inspect
 import numpy as np
@@ -477,7 +478,7 @@ def _format_parameters_table(parameter_names, fixed_width=None):
 
     if fixed_width is None:
         max_name_length = max([len(name) \
-                               for name in parameter_names + ['Name']])
+                               for name in parameter_names ]+[len('Name')])
         c1 = max_name_length + 1     # width of column 1
         c2 = (max_line_width - c1)   # width of column 2
     else:
@@ -521,8 +522,8 @@ def table_of_parameters(classname):
     opt_prm = getattr(classname, '_optional_parameters')
     for name in opt_prm:
         if not name in _parameters:
-            print 'Parameter "%s" used in class %s is not registered in _parameters.' % (name, classname.__name__)
-            print 'Do that before proceeding.'
+            print('Parameter "%s" used in class %s is not registered in _parameters.' % (name, classname.__name__))
+            print('Do that before proceeding.')
             sys.exit(1)
 
     s = """
@@ -696,10 +697,10 @@ def compile_string_functions(self, f, **kwargs):
                     else:
                         kwargs[func] = getattr(_callback,func)
             except:
-                raise ValueError, '''
+                raise ValueError('''
            F2py failed to compile input string (=\n%s)
            to be callable functions (%s).''' \
-                 % (string_to_compile, str_funcs.keys())
+                 % (string_to_compile, str_funcs.keys()))
         return f, kwargs
 
     def adjust_parameters(self):
@@ -859,13 +860,13 @@ def check_extra(self, **kwargs):
         for name, check_funcs, value in prm_type_list:
             try:
                 if not check_funcs(value):  # Return false
-                    raise ValueError,'''
+                    raise ValueError('''
         Improper value (=%s) for parameter %s.
-        Please check your input.''' % (str(value), name)
+        Please check your input.''' % (str(value), name))
             except:     # cannot run check_function smoothly
-                raise ValueError,'''
+                raise ValueError('''
         Cannot run check function for %s=%s.
-        Please check your input.''' % (name, str(value))
+        Please check your input.''' % (name, str(value)))
         return True
 
     def get(self, parameter_name=None, print_info=False):
@@ -883,20 +884,20 @@ def get(self, parameter_name=None, print_info=False):
             # user's functions. Instead, insert an entries that
             # reflect the name of user-supplied functions
             del all_args['f']
-            all_args['name of f'] = self.users_f.func_name
+            all_args['name of f'] = self.users_f.__name__
             if 'jac' in all_args:
                 del all_args['jac']
-                all_args['name of jac'] = self.users_jac.func_name
+                all_args['name of jac'] = self.users_jac.__name__
 
             if print_info:
-                print pprint.pformat(all_args)
+                print(pprint.pformat(all_args))
             return all_args
 
         else:
             if hasattr(self, parameter_name):
                 value = getattr(self, parameter_name)
                 if print_info:
-                    print "%s = %s" % (parameter_name, value)
+                    print("%s = %s" % (parameter_name, value))
                 return value
             else:
                 raise AttributeError('Parameter %s is not set' % parameter_name)
@@ -911,8 +912,8 @@ def get_parameter_info(self,print_info=False):
         is pretty printed, otherwise it is returned.
         '''
         if print_info:
-            print 'Legal parameters for class %s are:' % self.__class__.__name__
-            print pprint.pformat(self._parameters)
+            print('Legal parameters for class %s are:' % self.__class__.__name__)
+            print(pprint.pformat(self._parameters))
             return None
         else:
             return self._parameters
@@ -931,10 +932,10 @@ def _print_method(self, with_f, default):
         if with_f and hasattr(self, 'users_f'):
             if not hasattr(self.users_f, '__name__'):     # class instance?
                 f_name = self.users_f.__class__.__name__
-	    else:    # Ordinary functions
+            else:    # Ordinary functions
                 f_name = self.users_f.__name__
-	        if f_name == '<lambda>':   # lambda function
-	  	    f_name = 'lambda u, t: ...'
+            if f_name == '<lambda>':   # lambda function
+                f_name = 'lambda u, t: ...'
             args.append('f=%s' % f_name)
 
         # form all parameters
@@ -960,10 +961,10 @@ def _print_method(self, with_f, default):
         if hasattr(self, 'users_jac'):
             if hasattr(self.users_jac, '__name__'):     # plain function?
                 f_name = self.users_jac.__name__
-	        if f_name == '<lambda>':   # lambda function
-	  	    #f_name = 'lambda u, t: ...'
-	  	    f_name = 'lambda'
-	    else:                                       # class instance
+                if f_name == '<lambda>':   # lambda function
+                    #f_name = 'lambda u, t: ...'
+                    f_name = 'lambda'
+            else:                                       # class instance
                 f_name = self.users_jac.__class__.__name__
             args.append('jac=%s' % f_name)
 
@@ -1045,11 +1046,11 @@ def solve(self, time_points, terminate=None):
             self.u[n+1] = self.advance()   # new value
 
             if self.verbose > 2:
-                print '%s, step %d, t=%g' % \
-                      (self.__class__.__name__, n+1, self.t[n+1])
+                print('%s, step %d, t=%g' % \
+                      (self.__class__.__name__, n+1, self.t[n+1]))
             if terminate(self.u, self.t, n+1):
-                print self.__class__.__name__, \
-                      'terminated at t=%g' % self.t[n+1]
+                print(self.__class__.__name__, \
+                      'terminated at t=%g' % self.t[n+1])
                 if not self.disk_storage:
                     self.u, self.t = self.u[:n+2], self.t[:n+2]
                 # else: must keep original size of file, rest is 0
@@ -1084,7 +1085,7 @@ def initialize_for_solve(self):
             if len(self.t) < 2:
                 raise ValueError('time_points array %s must have at least two elements' % repr(self.t))
             if self.verbose > 0:
-                print 'Calling f(U0, %g) to determine data type' % self.t[0]
+                print('Calling f(U0, %g) to determine data type' % self.t[0])
             value = np.array(self.f(self.U0, self.t[0]))
         else:
             value = np.asarray(self.U0)
@@ -1161,23 +1162,20 @@ def validate_data(self):
             or (not np.asarray(
             # all items in self.t should be numbers
             [isinstance(t, (int,float)) for t in self.t]).all()):
-                raise TypeError, \
-                    'solve: time_points(=%s) is not a proper '\
-                    'sequence of real numbers' % str(self.t)
+                raise TypeError('solve: time_points(=%s) is not a proper '\
+                    'sequence of real numbers' % str(self.t))
 
         # self.t should be supplied in an asscending/descending order
         t_sorted = sorted(self.t, reverse=self.t[0] > self.t[-1])
         if list(self.t) != list(t_sorted):
-            raise ValueError, \
-                'time_points(=%s) is not provided in an ascending/descending'\
-                ' order!' % str(self.t)
+            raise ValueError('time_points(=%s) is not provided in an ascending/descending'\
+                ' order!' % str(self.t))
 
         # Test whether all required parameters are provided
         for arg in self._required_parameters:
             if not hasattr(self, arg):
-                raise ValueError,\
-                    '"%s" has to be input as required parameter(s) for '\
-                    'solver %s.' % (arg,self.__class__.__name__)
+                raise ValueError('"%s" has to be input as required parameter(s) for '\
+                    'solver %s.' % (arg,self.__class__.__name__))
         return True
 
     def switch_to(self, solver_target, print_info=False, **kwargs):
@@ -1220,10 +1218,10 @@ def switch_to(self, solver_target, print_info=False, **kwargs):
                 # Convert string to class name
                 solver_target = getattr(odespy, solver_target)
         except:
-            raise ValueError, error_message
+            raise ValueError(error_message)
 
         if not solver_target.__name__ in list_all_solvers():
-            raise ValueError, error_message
+            raise ValueError(error_message)
 
 
 
@@ -1266,9 +1264,9 @@ def switch_to(self, solver_target, print_info=False, **kwargs):
             diff_args = set(self.__dict__.keys()) - set(new.__dict__.keys()) \
                 - set(('u','t','n','dtype'))
             if diff_args:
-                print 'These attributes are neglected in %s: %s\n' \
-                    % (solver_target.__name__, str(diff_args)[5:-2])
-            print 'Switched to solver %s' % str(solver_target.__name__)
+                print('These attributes are neglected in %s: %s\n' \
+                    % (solver_target.__name__, str(diff_args)[5:-2]))
+            print('Switched to solver %s' % str(solver_target.__name__))
 
         return new
 
@@ -1334,9 +1332,9 @@ def check_conditional_parameters(self):
                         found = hasattr(self, arg)
                         arg_print = arg
                     if not found:
-                        raise ValueError,'''\
+                        raise ValueError('''\
         Error! Unsufficient input!
-        %s must be set when %s is %s!''' % (arg_print,name,value)
+        %s must be set when %s is %s!''' % (arg_print,name,value))
         return True
 
     def has_u_t_all(self):
@@ -1507,7 +1505,7 @@ def advance(self):
         dt = t[n+1] - t[n]
         K1 = dt*f(u[n], t[n])
         K2 = dt*f(u[n] + 0.5*K1, t[n] + 0.5*dt)
-	u_new = u[n] + K2
+        u_new = u[n] + K2
         return u_new
 
 
@@ -1587,7 +1585,7 @@ def initialize_for_solve(self):
     def validate_data(self):
         """Check that the time steps are constant."""
         if not self.constant_time_step():
-            print '%s must have constant time step' % self.__name__
+            print('%s must have constant time step' % self.__name__)
             return False
         else:
             return True
@@ -1636,7 +1634,7 @@ def initialize_for_solve(self):
 
     def validate_data(self):
         if not self.constant_time_step():
-            print '%s must have constant time step' % self.__name__
+            print('%s must have constant time step' % self.__name__)
             return False
         else:
             return True
@@ -1689,7 +1687,7 @@ def initialize_for_solve(self):
 
     def validate_data(self):
         if not self.constant_time_step():
-            print '%s must have constant time step' % self.__name__
+            print('%s must have constant time step' % self.__name__)
             return False
         else:
             return True
@@ -1743,7 +1741,7 @@ def initialize_for_solve(self):
 
     def validate_data(self):
         if not self.constant_time_step():
-            print '%s must have constant time step' % self.__name__
+            print('%s must have constant time step' % self.__name__)
             return False
         else:
             return True
@@ -1800,7 +1798,7 @@ def initialize_for_solve(self):
 
     def validate_data(self):
         if not self.constant_time_step():
-            print '%s must have constant time step' % self.__name__
+            print('%s must have constant time step' % self.__name__)
             return False
         else:
             return True
@@ -2013,7 +2011,7 @@ def initialize(self):
             import sympy
             self.sympy = sympy
         except ImportError:
-            raise ImportError,'sympy is not installed - needed for sympy_odefun'
+            raise ImportError('sympy is not installed - needed for sympy_odefun')
 
     def initialize_for_solve(self):
         # sympy.odefun requires f(t, u), not f(u, t, *args, **kwargs)
@@ -2119,8 +2117,8 @@ def advance(self):
         u_new = un + (t_new-tn)*f(un,tn)
         # control by number of intern steps and error tolerance
         if self.verbose > 1:
-            print '%s.advance w/%s: t=%g, n=%d: ' % \
-                  (self.__class__.__name__, self.nonlinear_solver, t_new, n+1),
+            print('%s.advance w/%s: t=%g, n=%d: ' % \
+                  (self.__class__.__name__, self.nonlinear_solver, t_new, n+1), end=' ')
 
         while i <= self.max_iter and error > self.eps_iter:
             if self.nonlinear_solver == 'Linear':
@@ -2140,11 +2138,11 @@ def advance(self):
             r = self.relaxation    # relaxation factor
             u_new = r*u_new_ + (1-r)*un
             if self.verbose > 2:
-                print '    u_%02d=%g (err=%g)' % (i, u_new, error)
+                print('    u_%02d=%g (err=%g)' % (i, u_new, error))
             i += 1
         self.num_iterations_total += i-1
         if self.verbose > 1:
-            print ' u=%g in %d iterations' % (u_new, i-1)
+            print(' u=%g in %d iterations' % (u_new, i-1))
         if error > self.eps_iter:
             raise ValueError('%s w/%s not converged:\n   difference in solution between last two iterations: %g > eps_iter=%g after %s iterations.' % (self.__class__.__name__, self.nonlinear_solver, error, self.eps_iter, self.max_iter))
         return u_new
@@ -2407,13 +2405,13 @@ def solve(self, time_points, terminate=None):
                 time_points = np.linspace(t[k-1], t[k], ntpoints+1)
                 dt = time_points[1] - time_points[0]
                 if self.verbose > 0:
-                    print 'dt=%g, ' % dt,
+                    print('dt=%g, ' % dt, end=' ')
                 if dt < self.min_step and self.verbose > 0:
-                    print 'too small %s < %s, ' % (dt, self.min_step)
-                    print
+                    print('too small %s < %s, ' % (dt, self.min_step))
+                    print()
                     break
                 if dt > self.max_step and self.verbose > 0:
-                    print 'too large %s > %s, ' % (dt, self.max_step)
+                    print('too large %s > %s, ' % (dt, self.max_step))
                     break
 
                 self.solver.set_initial_condition(self.u[-1])
@@ -2421,10 +2419,10 @@ def solve(self, time_points, terminate=None):
                 #R = self.residual(u_new[-2], u_new[-1], t_new[-2], t_new[-1])
                 R = self.residual(u_new[0], u_new[-1], t_new[0], t_new[-1])
                 if self.verbose > 0:
-                    print '\n%d time points in (t[%d], t[%d]) = (%.3g, %.3g), ' \
-                        % (ntpoints+1, k-1, k, t[k-1], t[k])
-                    print 'Residual=%.1E, tolerance=%.1E, calling %s' % \
-                        (R, self.atol, self.solver.__class__.__name__)
+                    print('\n%d time points in (t[%d], t[%d]) = (%.3g, %.3g), ' \
+                        % (ntpoints+1, k-1, k, t[k-1], t[k]))
+                    print('Residual=%.1E, tolerance=%.1E, calling %s' % \
+                        (R, self.atol, self.solver.__class__.__name__))
                 # reintegrate with time_step = dt/2
                 ntpoints *= 2
             self.u.extend(u_new[1:])
@@ -2511,8 +2509,8 @@ def middle(x, y, z):
         self.h = first_step
 
         if self.verbose:
-            print '\nadvance solution in [%g, %g], ' % (t, t_np1),
-            print 'starting with h=%g' % self.h
+            print('\nadvance solution in [%g, %g], ' % (t, t_np1), end=' ')
+            print('starting with h=%g' % self.h)
 
         while t < t_np1:
 
@@ -2521,7 +2519,7 @@ def middle(x, y, z):
                 u_new, error = self.advance_intermediate(u, t, self.h)
 
                 if self.verbose:
-                    print 'u(t=%g): %g, ' % (t + self.h, u_new),
+                    print('u(t=%g): %g, ' % (t + self.h, u_new), end=' ')
 
                 # Is u_new sufficiently accurate?
                 u_new_norm = np.linalg.norm(u_new)
@@ -2537,24 +2535,24 @@ def middle(x, y, z):
 
                 if self.verbose:
                     if sufficiently_accurate:
-                        print 'accepted, ',
+                        print('accepted, ', end=' ')
                     else:
-                        print 'rejected, ',
-                    print ' err=%g (tol=%.1E), ' % (error, tol),
+                        print('rejected, ', end=' ')
+                    print(' err=%g (tol=%.1E), ' % (error, tol), end=' ')
                     if hasattr(self, 'u_exact'):
-                        print 'exact-err: %g, ' % \
-                              (np.abs(np.asarray(self.u_exact(t))-u)),
+                        print('exact-err: %g, ' % \
+                              (np.abs(np.asarray(self.u_exact(t))-u)), end=' ')
 
                 # Adjust time step
                 if error > 1E-14:
                     h_new = self.h*(tol/error)**(1./self.order)
-                    print 'changing time step', h_new, min_step, max_step
+                    print('changing time step', h_new, min_step, max_step)
                     self.h = middle(min_step, h_new, max_step)
 
                 self.h = min(self.h, t_np1-t)  # fit last step
 
                 if self.verbose:
-                    print 'new h=%g' % self.h
+                    print('new h=%g' % self.h)
         return u
 
 
@@ -2727,7 +2725,7 @@ def advance(self):
         u, f, n, t = self.u, self.f, self.n, self.t
         u_new = self.integrator.integrate(t[n+1])
         if not self.integrator.successful():
-            print 'Warning: %s call to ode.method.integrate in scipy.integrate was not successful' % self.__class__.__name__
+            print('Warning: %s call to ode.method.integrate in scipy.integrate was not successful' % self.__class__.__name__)
         if len(u_new) == 1:
             return u_new[0]
         else:
@@ -2908,7 +2906,7 @@ class odelab(Adaptive):
     solvers = 'ExplicitEuler ExplicitTrapezoidal ImplicitEuler RungeKutta34 RungeKutta4 SymplecticEuler ImplicitEuler LDIRK343 LobattoIIIA LobattoIIIB LobattoIIIC LobattoIIICs LobattoIIID MultiRKDAE RKDAE RadauIIA RungeKutta Spark Spark2 Spark3 Spark4 AdamsBashforth AdamsBashforth1 AdamsBashforth2 AdamsBashforth2e Butcher Butcher1 Butcher3 ExplicitGeneralLinear GeneralLinear Heun Kutta Kutta38 Kutta4 McLachlan NonHolonomicEnergy NonHolonomicEnergy0 NonHolonomicEnergyEMP NonHolonomicLeapFrog SymplecticEuler ABLawson ABLawson2 ABLawson3 ABLawson4 ABNorset4 Exponential GenLawson45 HochOst4 Lawson4 LawsonEuler Phi Polynomial RKMK4T'.split()
     DAE_solvers = 'SymplicticEuler MultiRKDAE RKDAE Spark Spark2 Spark3 Spark4'.split()
     not_valid = 'SymplecticEuler LDIRK343 LobattoIIIA LobattoIIIB LobattoIIIC LobattoIIICs LobattoIIID MultiRKDAE RKDAE RadauIIA RungeKutta Spark Spark2 Spark3 Spark4 AdamsBashforth AdamsBashforth2 AdamsBashforth2e Butcher Butcher1 Butcher3 ExplicitGeneralLinear GeneralLinear Kutta McLachlan NonHolonomicEnergy NonHolonomicEnergy0 NonHolonomicEnergyEMP NonHolonomicLeapFrog SymplecticEuler ABLawson ABLawson2 ABLawson3 ABLawson4 ABNorset4 Exponential GenLawson45 HochOst4 Lawson4 LawsonEuler Phi Polynomial RKMK4T'.split()
-    solvers = [solver for solver in solvers if not solver in not_valid]
+    solvers = list(set(solvers) - set(not_valid))
 
     __doc__ = __doc__ % (str(solvers)[1:-1])
 
@@ -2943,7 +2941,7 @@ def initialize(self):
             self.odelab_schemes = schemes
 
         except ImportError:
-            raise ImportError,'The odelab software is not installed - needed for class odelab'
+            raise ImportError('The odelab software is not installed - needed for class odelab')
 
     def initialize_for_solve(self):
         # odelab requires f(t, u), not f(u, t, *args, **kwargs)
@@ -2979,7 +2977,7 @@ def solve(self, time_points, terminate=None):
         since odelab. is such a solve method.
         """
         if terminate is not None:
-            print 'Warning: odelab.solve ignores the terminate function!'
+            print('Warning: odelab.solve ignores the terminate function!')
         self.t = np.asarray(time_points)
 
         try:
@@ -2992,7 +2990,7 @@ def solve(self, time_points, terminate=None):
             self.solver.initialize(self.U0)
             # Solve problem in odelab
             self.solver.run(self.t[-1])
-        except Exception, e:
+        except Exception as e:
             msg = 'odelab scheme %s did not work:\n%s' % \
                   (self.odelab_solver, e)
             #print msg
diff --git a/odespy/tests/joblib_ex.py b/odespy/tests/joblib_ex.py
index 3256755..8916163 100644
--- a/odespy/tests/joblib_ex.py
+++ b/odespy/tests/joblib_ex.py
@@ -28,13 +28,14 @@
 Is this the right way to use joblib? Or is the same functionality available
 in better ways?
 """
+from __future__ import print_function
 
 from joblib import Memory
 memory = Memory(cachedir='tmp', verbose=0)
 
 @memory.cache(ignore=['data'])
 def retrieve(name, data=None):
-    print 'joblib save of', name
+    print('joblib save of', name)
     return data
 
 save = retrieve
@@ -70,7 +71,7 @@ def compute(data):
     create_big_data(1000)  # create pieces and store in joblib
     data = retrieve_big_data()
     data = compute(data)
-    print data['arrays']['b'][-12:]
+    print(data['arrays']['b'][-12:])
 
 # Try to do different things with b and see if retrieve is ever
 # called.
diff --git a/odespy/tests/test_ThetaRule.py b/odespy/tests/test_ThetaRule.py
index 6a9100d..bcb3661 100644
--- a/odespy/tests/test_ThetaRule.py
+++ b/odespy/tests/test_ThetaRule.py
@@ -1,3 +1,4 @@
+from __future__ import print_function
 from numpy.testing import assert_array_almost_equal, TestCase, \
     run_module_suite, assert_almost_equal
 
@@ -64,9 +65,9 @@ def myjac(u, t):
                              rtol=10*prms0['eps_iter'],
                              atol=prms0['eps_iter'])
             # str(solver) prints all parameters *different* from the defaults
-            print str(solver), 'works' if ok else 'does not work'
+            print(str(solver), 'works' if ok else 'does not work')
             if not ok:
-                print 'max deviation:', abs(u_d-u).max()
+                print('max deviation:', abs(u_d-u).max())
             assert ok
 
 if __name__ == "__main__":
diff --git a/odespy/tests/test_basics.py b/odespy/tests/test_basics.py
index 7b16fd4..4f3ac23 100644
--- a/odespy/tests/test_basics.py
+++ b/odespy/tests/test_basics.py
@@ -12,6 +12,7 @@
 future simulations with the recorded difference (check that these are
 within machine precision, here taken as 1E-14).
 """
+from __future__ import print_function
 # Author: Liwei Wang and Hans Petter Langtangen
 
 import odespy
@@ -84,8 +85,8 @@
 
 Sine = dict(
     help="Sine, u'' = -u, u = sin(t)",
-    f=lambda (u00,u11),t : [u11, -u00],
-    jac=lambda (u00,u11),t: [[0., 1.], [- 1., 0.]],
+    f=lambda u,t : [u[1], -u[0]],
+    jac=lambda u,t: [[0., 1.], [- 1., 0.]],
     time_points=np.linspace(0., 1, 10),
     # These solvers are not suitable for this ODE problem
     exceptions=['Lsodi', 'Lsodes', 'Lsodis', 'Lsoibt', 'MyRungeKutta',
@@ -142,8 +143,8 @@
 
 VanDerPol = dict(
     help="Van der Pol oscillator problem:  u'' = 3*(1 - u**2)*u' - u",
-    f=lambda (u00,u11),t: [u11, 3.*(1 - u00**2)*u11 - u00],
-    jac=lambda (u00,u11),t: [[0., 1.], [-6.*u00*u11 - 1., 3.*(1. - u00**2)]],
+    f=lambda u,t: [u[1], 3.*(1 - u[0]**2)*u[1] - u[0]],
+    jac=lambda u,t: [[0., 1.], [-6.*u[0]*u[1] - 1., 3.*(1. - u[0]**2)]],
     u0=[2., 0.],
     # These solvers are not suitable for this ODE problem
     exceptions=['RKC', 'Lsodes', 'Leapfrog', 'Lsodi', 'Lsodis', 'Lsoibt',
@@ -244,11 +245,11 @@ def _run_test_problems(problem):
                      'f_kwargs')
     kwargs = {key: problem[key] for key in problem
               if key not in special_names}
-    print problem.get('help', '')
+    print(problem.get('help', ''))
     for method in methods:
 
         solver = eval('odespy.%s' % method)(problem['f'], **kwargs)
-        print 'Testing %s' % method,
+        print('Testing %s' % method, end=' ')
         solver.set_initial_condition(problem['u0'])
         u, t = solver.solve(problem['time_points'])
 
@@ -277,10 +278,10 @@ def _run_test_problems(problem):
             if method in problem['exact_diff']:
                 nt.assert_almost_equal(diff, problem['exact_diff'][method],
                                        delta=1E-14)
-                print '...ok'
+                print('...ok')
             else:
                 pass
-                print ' no exact diff available for comparison'
+                print(' no exact diff available for comparison')
         if 'exact_final_diff' in problem:
             u_final = np.asarray(u[-1])
             exact_final = np.asarray(problem['exact_final'])
@@ -288,9 +289,9 @@ def _run_test_problems(problem):
             if method in problem['exact_final_diff']:
                 nt.assert_almost_equal(diff, problem['exact_final_diff'][method],
                                        delta=1E-14)
-                print '...ok'
+                print('...ok')
             else:
-                print ' no exact final diff available for comparison'
+                print(' no exact final diff available for comparison')
 
 def test_exponentinal():
     _run_test_problems(Exponential)
@@ -306,7 +307,7 @@ def test_complex():
 
 def test_EulerCromer_1dof():
     """Plain u'' + u = 0."""
-    solver = odespy.EulerCromer(lambda (v, x), t: [-x, v])
+    solver = odespy.EulerCromer(lambda vx, t: [-vx[1], vx[0]])
     #solver = odespy.EulerCromer(lambda u, t: [-u[1], u[0]])
     solver.set_initial_condition([0, 1])
     P = 2*np.pi
@@ -321,7 +322,7 @@ def test_EulerCromer_1dof():
     x_exact = lambda t: np.cos(t)
     x_e = x_exact(t)
     #plot(t, x, t, x_e, legend=('EC', 'exact'))
-    print 'Testing EulerCromer',
+    print('Testing EulerCromer', end=' ')
     # Test difference in final value
     if N == 60:
         diff_exact = 0.0014700112828
@@ -329,7 +330,7 @@ def test_EulerCromer_1dof():
         tol = 1E-14
         assert abs(diff_exact - diff_EC) < tol, \
                'diff_exact=%g, diff_EulerCromer=%g' % (diff_exact, diff_EC)
-        print '...ok'
+        print('...ok')
 
 def test_terminate():
     problem = Exponential
@@ -338,9 +339,9 @@ def test_terminate():
     u, t = solver.solve(problem['time_points'], terminate=problem['terminate'])
     exact_diff = 0.0321206486802586
     diff = abs(problem['stop_value']-u[-1])
-    print 'Testing RKFehlberg with terminate function',
+    print('Testing RKFehlberg with terminate function', end=' ')
     nt.assert_almost_equal(diff, exact_diff, delta=1E-14)
-    print '...ok'
+    print('...ok')
 
 def test_switch_to():
     problem = Exponential
@@ -351,9 +352,9 @@ def test_switch_to():
     u2, t2 = solver_new.solve(problem['time_points'])
     diff = np.abs(u - u2).max()
     exact_diff = 0.0000015284633990
-    print 'Testing switch_to from RKFehlberg to Lsode',
+    print('Testing switch_to from RKFehlberg to Lsode', end=' ')
     nt.assert_almost_equal(diff, exact_diff, delta=1E-14)
-    print '...ok'
+    print('...ok')
 
 
 if __name__ == '__main__':