diff --git a/.travis.yml b/.travis.yml
index 65cdc8550..a28f57419 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -76,5 +76,18 @@ script:
   - 'if [ $SLYCOT != "" ]; then python -c "import slycot"; fi'
   - coverage run setup.py test
 
+  # only run examples if Slycot is install
+  # set PYTHONPATH for examples
+  # pmw needed for examples/tfvis.py
+  # future is needed for Python 2, also for examples/tfvis.py
+
+  - if [[ "$SLYCOT" != "" ]]; then
+      export PYTHONPATH=$PWD;
+      conda install -c conda-forge pmw future;
+      cd examples; bash run_examples.sh; cd ..;
+    fi
+
+# arbitrary change to try to trigger travis build
+
 after_success:
   - coveralls
diff --git a/control/phaseplot.py b/control/phaseplot.py
index a434985ae..10fbec640 100644
--- a/control/phaseplot.py
+++ b/control/phaseplot.py
@@ -39,12 +39,20 @@
 
 import numpy as np
 import matplotlib.pyplot as mpl
-from matplotlib.mlab import frange, find
+
 from scipy.integrate import odeint
 from .exception import ControlNotImplemented
 
 __all__ = ['phase_plot', 'box_grid']
 
+
+def _find(condition):
+    """Returns indices where ravel(a) is true.
+    Private implementation of deprecated matplotlib.mlab.find
+    """
+    return np.nonzero(np.ravel(condition))[0]
+
+
 def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None,
               lingrid=None, lintime=None, logtime=None, timepts=None,
               parms=(), verbose=True):
@@ -70,11 +78,11 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None,
         dxdt = F(x, t) that accepts a state x of dimension 2 and
         returns a derivative dx/dt of dimension 2.
 
-    X, Y: ndarray, optional
-        Two 1-D arrays representing x and y coordinates of a grid.
-        These arguments are passed to meshgrid and generate the lists
-        of points at which the vector field is plotted.  If absent (or
-        None), the vector field is not plotted.
+    X, Y: 3-element sequences, optional, as [start, stop, npts]
+        Two 3-element sequences specifying x and y coordinates of a
+        grid.  These arguments are passed to linspace and meshgrid to
+        generate the points at which the vector field is plotted.  If
+        absent (or None), the vector field is not plotted.
 
     scale: float, optional
         Scale size of arrows; default = 1
@@ -145,8 +153,10 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None,
     #! TODO: Add sanity checks
     elif (X is not None and Y is not None):
         (x1, x2) = np.meshgrid(
-            frange(X[0], X[1], float(X[1]-X[0])/X[2]),
-            frange(Y[0], Y[1], float(Y[1]-Y[0])/Y[2]));
+            np.linspace(X[0], X[1], X[2]),
+            np.linspace(Y[0], Y[1], Y[2]))
+        Narrows = len(x1)
+
     else:
         # If we weren't given any grid points, don't plot arrows
         Narrows = 0;
@@ -234,12 +244,12 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None,
                 elif (logtimeFlag):
                     # Use an exponential time vector
                     # MATLAB: tind = find(time < (j-k) / lambda, 1, 'last');
-                    tarr = find(time < (j-k) / timefactor);
+                    tarr = _find(time < (j-k) / timefactor);
                     tind = tarr[-1] if len(tarr) else 0;
                 elif (timeptsFlag):
                     # Use specified time points
                     # MATLAB: tind = find(time < Y[j], 1, 'last');
-                    tarr = find(time < timepts[j]);
+                    tarr = _find(time < timepts[j]);
                     tind = tarr[-1] if len(tarr) else 0;
 
                 # For tailless arrows, skip the first point
@@ -295,8 +305,8 @@ def box_grid(xlimp, ylimp):
     box defined by the corners [xmin ymin] and [xmax ymax].
     """
 
-    sx10 = frange(xlimp[0], xlimp[1], float(xlimp[1]-xlimp[0])/xlimp[2])
-    sy10 = frange(ylimp[0], ylimp[1], float(ylimp[1]-ylimp[0])/ylimp[2])
+    sx10 = np.linspace(xlimp[0], xlimp[1], xlimp[2])
+    sy10 = np.linspace(ylimp[0], ylimp[1], ylimp[2])
 
     sx1 = np.hstack((0, sx10, 0*sy10+sx10[0], sx10, 0*sy10+sx10[-1]))
     sx2 = np.hstack((0, 0*sx10+sy10[0], sy10, 0*sx10+sy10[-1], sy10))
diff --git a/examples/genswitch.py b/examples/genswitch.py
index 7e36fc357..11f79a36c 100644
--- a/examples/genswitch.py
+++ b/examples/genswitch.py
@@ -5,10 +5,11 @@
 # of a genetic switch.  Plots time traces and a phase portrait using
 # the python-control library.
 
+import os
+
 import numpy as np
 import matplotlib.pyplot as mpl
 from scipy.integrate import odeint
-from matplotlib.mlab import frange
 from control import phase_plot, box_grid
 
 # Simple model of a genetic switch
@@ -34,7 +35,7 @@ def genswitch(y, t, mu=4, n=2):
 sol2 = odeint(genswitch, sol1[-1,:] + [2, -2], tim2)
 
 # First plot out the curves that define the equilibria
-u = frange(0, 4.5, 0.1)
+u = np.linspace(0, 4.5, 46)
 f = np.divide(mu, (1 + u**n))   # mu / (1 + u^n), elementwise
 
 mpl.figure(1); mpl.clf();
@@ -51,7 +52,6 @@ def genswitch(y, t, mu=4, n=2):
 mpl.figure(3); mpl.clf(); # subplot(221);
 mpl.plot(tim1, sol1[:,0], 'b-', tim1, sol1[:,1], 'g--');
 # set(pl, 'LineWidth', AM_data_linewidth);
-mpl.hold(True);
 mpl.plot([tim1[-1], tim1[-1]+1], 
          [sol1[-1,0], sol2[0,1]], 'ko:', 
          [tim1[-1], tim1[-1]+1], [sol1[-1,1], sol2[0,0]], 'ko:');
@@ -72,10 +72,11 @@ def genswitch(y, t, mu=4, n=2):
           timepts = [0.2, 0.6, 1.2])
 
 # Add the stable equilibrium points
-mpl.hold(True);
 mpl.plot(eqpt[0], eqpt[1], 'k.', eqpt[1], eqpt[0], 'k.', 
   eqpt[2], eqpt[2], 'k.')       # 'MarkerSize', AM_data_markersize*3);
 
 mpl.xlabel('Protein A [scaled]');
 mpl.ylabel('Protein B [scaled]');       # 'Rotation', 90);
 
+if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
+    mpl.show()
diff --git a/examples/phaseplots.py b/examples/phaseplots.py
index 3d74dca7d..0c86522de 100644
--- a/examples/phaseplots.py
+++ b/examples/phaseplots.py
@@ -4,6 +4,8 @@
 # This file contains examples of phase portraits pulled from "Feedback
 # Systems" by Astrom and Murray (Princeton University Press, 2008).
 
+import os
+
 import numpy as np
 import matplotlib.pyplot as mpl
 from control.phaseplot import phase_plot
@@ -23,7 +25,7 @@ def invpend_ode(x, t, m=1., l=1., b=0.2, g=1):
 # Set up the figure the way we want it to look
 mpl.figure(); mpl.clf(); 
 mpl.axis([-2*pi, 2*pi, -2.1, 2.1]);
-mpl.title('Inverted pendlum')
+mpl.title('Inverted pendulum')
 
 # Outer trajectories
 phase_plot(invpend_ode,
@@ -35,9 +37,7 @@ def invpend_ode(x, t, m=1., l=1., b=0.2, g=1):
     logtime = (3, 0.7) )
 
 # Separatrices
-mpl.hold(True);
 phase_plot(invpend_ode, X0 = [[-2.3056, 2.1], [2.3056, -2.1]], T=6, lingrid=0)
-mpl.show();
 
 #
 # Systems of ODEs: damped oscillator example (simulation + phase portrait)
@@ -49,9 +49,10 @@ def oscillator_ode(x, t, m=1., b=1, k=1):
 # Generate a vector plot for the damped oscillator
 mpl.figure(); mpl.clf();
 phase_plot(oscillator_ode, [-1, 1, 10], [-1, 1, 10], 0.15);
-mpl.hold(True); mpl.plot([0], [0], '.');
+#mpl.plot([0], [0], '.');
 # a=gca; set(a,'FontSize',20); set(a,'DataAspectRatio',[1,1,1]);
-mpl.xlabel('x1'); mpl.ylabel('x2');
+mpl.xlabel('$x_1$'); mpl.ylabel('$x_2$');
+mpl.title('Damped oscillator, vector field')
 
 # Generate a phase plot for the damped oscillator
 mpl.figure(); mpl.clf(); 
@@ -61,11 +62,10 @@ def oscillator_ode(x, t, m=1., b=1, k=1):
     [-1, 1], [-0.3, 1], [0, 1], [0.25, 1], [0.5, 1], [0.75, 1], [1, 1],
     [1, -1], [0.3, -1], [0, -1], [-0.25, -1], [-0.5, -1], [-0.75, -1], [-1, -1]
   ], T = np.linspace(0, 8, 80), timepts = [0.25, 0.8, 2, 3])
-mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
+mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
 # set(gca,'DataAspectRatio',[1,1,1]);
-mpl.xlabel('x1'); mpl.ylabel('x2');
-
-mpl.show()
+mpl.xlabel('$x_1$'); mpl.ylabel('$x_2$');
+mpl.title('Damped oscillator, vector field and stream lines')
 
 #
 # Stability definitions
@@ -88,9 +88,10 @@ def saddle_ode(x, t):
     [-1.3,-1]
   ], T = np.linspace(0, 10, 100), 
   timepts = [0.3, 1, 2, 3], parms = (m, b, k));
-mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
+mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
 # set(gca,'FontSize', 16); 
-mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
+mpl.xlabel('$x_1$'); mpl.ylabel('$x_2$');
+mpl.title('Asymptotically stable point')
 
 # Saddle
 mpl.figure(); mpl.clf();
@@ -103,9 +104,10 @@ def saddle_ode(x, t):
     [0.95, 1], [0.9, 1], [0.8, 1], [0.6, 1], [0.4, 1], [0.2, 1],
     [-0.5, -0.45], [-0.45, -0.5], [0.5, 0.45], [0.45, 0.5],
     [-0.04, 0.04], [0.04, -0.04] ], T = np.linspace(0, 2, 20));
-mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
+mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
 # set(gca,'FontSize', 16); 
-mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
+mpl.xlabel('$x_1$'); mpl.ylabel('$x_2$');
+mpl.title('Saddle point')
 
 # Stable isL
 m = 1; b = 0; k = 1;			# zero damping
@@ -115,8 +117,10 @@ def saddle_ode(x, t):
   [pi/6, pi/3, pi/2, 2*pi/3, 5*pi/6, pi, 7*pi/6, 4*pi/3, 9*pi/6, 5*pi/3, 11*pi/6, 2*pi], 
   X0 = [ [0.2,0], [0.4,0], [0.6,0], [0.8,0], [1,0], [1.2,0], [1.4,0] ],
   T = np.linspace(0, 20, 200), parms = (m, b, k));
-mpl.hold(True); mpl.plot([0], [0], 'k.') # 'MarkerSize', AM_data_markersize*3);
+mpl.plot([0], [0], 'k.') # 'MarkerSize', AM_data_markersize*3);
 # set(gca,'FontSize', 16); 
-mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
+mpl.xlabel('$x_1$'); mpl.ylabel('$x_2$');
+mpl.title('Undamped system\nLyapunov stable, not asympt. stable')
 
-mpl.show()
+if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
+    mpl.show()
diff --git a/examples/pvtol-lqr.py b/examples/pvtol-lqr.py
index d50584163..8412dc2ff 100644
--- a/examples/pvtol-lqr.py
+++ b/examples/pvtol-lqr.py
@@ -3,10 +3,12 @@
 #
 # This file works through an LQR based design problem, using the
 # planar vertical takeoff and landing (PVTOL) aircraft example from
-# Astrom and Mruray, Chapter 5.  It is intended to demonstrate the
+# Astrom and Murray, Chapter 5.  It is intended to demonstrate the
 # basic functionality of the python-control package.
 #
 
+import os
+
 from numpy import *             # Grab all of the NumPy functions
 from matplotlib.pyplot import * # Grab MATLAB plotting functions
 from control.matlab import *    # MATLAB-like functions
@@ -118,9 +120,9 @@
 (Yy, Ty) = step(H1ay, T=linspace(0,10,100));
 
 subplot(221); title("Identity weights")
-# plot(T, Y[:,1, 1], '-', T, Y[:,2, 2], '--'); hold(True);
-plot(Tx.T, Yx.T, '-', Ty.T, Yy.T, '--'); hold(True);
-plot([0, 10], [1, 1], 'k-'); hold(True);
+# plot(T, Y[:,1, 1], '-', T, Y[:,2, 2], '--');
+plot(Tx.T, Yx.T, '-', Ty.T, Yy.T, '--'); 
+plot([0, 10], [1, 1], 'k-'); 
 
 axis([0, 10, -0.1, 1.4]); 
 ylabel('position');
@@ -141,10 +143,10 @@
 [Y3, T3] = step(H1cx, T=linspace(0,10,100));
 
 subplot(222); title("Effect of input weights")
-plot(T1.T, Y1.T, 'b-'); hold(True);
-plot(T2.T, Y2.T, 'b-'); hold(True);
-plot(T3.T, Y3.T, 'b-'); hold(True);
-plot([0 ,10], [1, 1], 'k-'); hold(True);
+plot(T1.T, Y1.T, 'b-');
+plot(T2.T, Y2.T, 'b-');
+plot(T3.T, Y3.T, 'b-');
+plot([0 ,10], [1, 1], 'k-');
 
 axis([0, 10, -0.1, 1.4]); 
 
@@ -190,4 +192,6 @@
 xlabel('time'); 
 legend(('x', 'y'), loc='lower right');
 
-show()
+if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
+    show()
+
diff --git a/examples/pvtol-nested-ss.py b/examples/pvtol-nested-ss.py
index d78e97981..24e173bc8 100644
--- a/examples/pvtol-nested-ss.py
+++ b/examples/pvtol-nested-ss.py
@@ -98,7 +98,7 @@
 bode(L, logspace(-4, 3));
 
 # Add crossover line
-subplot(211); hold(True);
+subplot(211);
 loglog([1e-4, 1e3], [1, 1], 'k-')
 
 # Replot phase starting at -90 degrees
@@ -107,7 +107,6 @@
 
 subplot(212);
 semilogx([1e-4, 1e3], [-180, -180], 'k-')
-hold(True);
 semilogx(w, np.squeeze(phase), 'b-')
 axis([1e-4, 1e3, -360, 0]);
 xlabel('Frequency [deg]'); ylabel('Phase [deg]');
@@ -118,7 +117,7 @@
 # Nyquist plot for complete design
 #
 figure(7); clf;
-axis([-700, 5300, -3000, 3000]); hold(True);
+axis([-700, 5300, -3000, 3000]);
 nyquist(L, (0.0001, 1000));
 axis([-700, 5300, -3000, 3000]);
 
@@ -127,7 +126,7 @@
 
 # Expanded region  
 figure(8); clf; subplot(231); 
-axis([-10, 5, -20, 20]); hold(True);
+axis([-10, 5, -20, 20]);
 nyquist(L);
 axis([-10, 5, -20, 20]);
 
@@ -145,7 +144,7 @@
 
 figure(9); 
 (Yvec, Tvec) = step(T, linspace(1, 20));
-plot(Tvec.T, Yvec.T); hold(True);
+plot(Tvec.T, Yvec.T);
 
 (Yvec, Tvec) = step(Co*S, linspace(1, 20));
 plot(Tvec.T, Yvec.T);
diff --git a/examples/robust_mimo.py b/examples/robust_mimo.py
index 5270579d8..15abc9263 100644
--- a/examples/robust_mimo.py
+++ b/examples/robust_mimo.py
@@ -4,6 +4,8 @@
 and Postlethwaite, 1st Edition.
 """
 
+import os
+
 import numpy as np
 import matplotlib.pyplot as plt
 
@@ -68,7 +70,7 @@ def analysis():
     plt.ylabel('output')
     plt.ylim([-1.1,2.1])
     plt.legend()
-    plt.title('o/l response to input [1,0]')
+    plt.title('o/l response\nto input [1,0]')
 
     plt.subplot(1,3,2)
     plt.plot(t,yu2[0],label='$y_1$')
@@ -77,7 +79,7 @@ def analysis():
     plt.ylabel('output')
     plt.ylim([-1.1,2.1])
     plt.legend()
-    plt.title('o/l response to input [0,1]')
+    plt.title('o/l response\nto input [0,1]')
 
     plt.subplot(1,3,3)
     plt.plot(t,yuz[0],label='$y_1$')
@@ -86,7 +88,7 @@ def analysis():
     plt.ylabel('output')
     plt.ylim([-1.1,2.1])
     plt.legend()
-    plt.title('o/l response to input [1,-1]')
+    plt.title('o/l response\nto input [1,-1]')
 
 
 def synth(wb1,wb2):
@@ -103,7 +105,7 @@ def synth(wb1,wb2):
     wp2 = ss(weighting(wb=wb2, m=1.5, a=1e-4))
     wp = wp1.append(wp2)
     k,_,info = mixsyn(g,wp,wu)
-    return k, info.gamma
+    return k, info[0]
 
 
 def step_opposite(g,t):
@@ -177,4 +179,6 @@ def design():
 
 analysis()
 design()
-plt.show()
+if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
+    plt.show()
+
diff --git a/examples/robust_siso.py b/examples/robust_siso.py
index 39945919c..6bc4ab327 100644
--- a/examples/robust_siso.py
+++ b/examples/robust_siso.py
@@ -4,6 +4,8 @@
 and Postlethwaite, 1st Edition.
 """
 
+import os
+
 import numpy as np
 import matplotlib.pyplot as plt
 
@@ -99,4 +101,5 @@
 plt.legend()
 plt.title('Disturbance response')
 
-plt.show()
+if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
+    plt.show()
diff --git a/examples/rss-balred.py b/examples/rss-balred.py
index 7cadd043c..86e499a80 100755
--- a/examples/rss-balred.py
+++ b/examples/rss-balred.py
@@ -1,5 +1,7 @@
 #!/usr/bin/env python
 
+import os
+
 import numpy as np
 import control.modelsimp as msimp
 import control.matlab as mt
@@ -29,7 +31,6 @@
 y, t = mt.step(fsys)
 yr, tr = mt.step(rsys)
 plt.plot(t.T, y.T)
-plt.hold(True)
 plt.plot(tr.T, yr.T)
 
 # Repeat balanced reduction, now with 100-dimensional random state space
@@ -42,3 +43,7 @@
 yrandr, trandr = mt.impulse(rsysrand)
 plt.plot(trand.T, yrand.T, trandr.T, yrandr.T) 
 
+
+if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
+    plt.show()
+
diff --git a/examples/run_examples.sh b/examples/run_examples.sh
new file mode 100755
index 000000000..6f04fe12c
--- /dev/null
+++ b/examples/run_examples.sh
@@ -0,0 +1,27 @@
+#!/bin/bash
+# Run the examples.  Doesn't run Jupyter notebooks.
+
+# The examples must cooperate: they must have no operations that wait
+# on user action, like matplotlib.pyplot.show, or starting a GUI event
+# loop.  Environment variable PYCONTROL_TEST_EXAMPLES is set when the
+# examples are being tested; existence of this variable should be used
+# to prevent such user-action-waiting operations.
+
+export PYCONTROL_TEST_EXAMPLES=1
+
+example_errors=""
+
+for example in *.py; do
+    echo "Running ${example}"
+    if ! python ${example}; then
+        example_errors="${example_errors} ${example}"
+    fi
+done
+    
+if [ -n "${example_errors}" ]; then
+    echo These examples had errors:
+    echo "${example_errors}"
+    exit 1
+else
+    echo All examples ran successfully
+fi
diff --git a/examples/slicot-test.py b/examples/slicot-test.py
index f14e2995c..bf8af4206 100644
--- a/examples/slicot-test.py
+++ b/examples/slicot-test.py
@@ -1,6 +1,8 @@
 # Simple test function to test out SLICOT interface
 # RMM, 28 May 09
 
+from __future__ import print_function
+
 import numpy as np              # Numerical library
 from scipy import *             # Load the scipy functions
 from control.matlab import *    # Load the controls systems library
@@ -19,8 +21,8 @@
 # Eigenvalue placement
 #from slycot import sb01bd
 K = place(A, B, [-3, -2, -1])
-print "Pole place: K = ", K
-print "Pole place: eigs = ", np.linalg.eig(A - B * K)[0]
+print("Pole place: K = ", K)
+print("Pole place: eigs = ", np.linalg.eig(A - B * K)[0])
 
 # from slycot import ab01md
 # Ac, Bc, ncont, info = ab01md('I', A, B)
diff --git a/examples/test-statefbk.py b/examples/test-statefbk.py
index b7467bc7e..d9236d949 100644
--- a/examples/test-statefbk.py
+++ b/examples/test-statefbk.py
@@ -1,6 +1,8 @@
 # test-statefbk.py - Unit tests for state feedback code
 # RMM, 6 Sep 2010
 
+from __future__ import print_function
+
 import numpy as np              # Numerical library
 from scipy import *             # Load the scipy functions
 from control.matlab import *    # Load the controls systems library
@@ -18,10 +20,10 @@
 
 # Controllability
 Wc = ctrb(A, B)
-print "Wc = ", Wc
+print("Wc = ", Wc)
 
 # Observability
 Wo = obsv(A, C)
-print "Wo = ", Wo
+print("Wo = ", Wo)
 
 
diff --git a/examples/tfvis.py b/examples/tfvis.py
index 7929f7023..056fd62eb 100644
--- a/examples/tfvis.py
+++ b/examples/tfvis.py
@@ -1,4 +1,9 @@
 #!/usr/bin/python
+# needs pmw (in pypi, conda-forge)
+# For Python 2, needs future (in conda pypi and "default")
+
+from __future__ import print_function
+
 """ Simple GUI application for visualizing how the poles/zeros of the transfer
 function effects the bode, nyquist and step response of a SISO system """
 
@@ -38,7 +43,7 @@
 """
 
 import control.matlab
-import Tkinter
+import tkinter
 import sys
 import Pmw
 import matplotlib.pyplot as plt
@@ -101,7 +106,7 @@ def __init__(self, parent):
             pass
 
         widgets = (self.numerator_widget, self.denominator_widget)
-        for i in xrange(len(widgets)):
+        for i in range(len(widgets)):
             widgets[i].grid(row=i+1, column=0, padx=20, pady=3)
         Pmw.alignlabels(widgets)
 
@@ -156,16 +161,16 @@ def __init__(self, parent):
         self.zeros = []
         self.poles = []
 
-        self.topframe = Tkinter.Frame(self.master)
+        self.topframe = tkinter.Frame(self.master)
         self.topframe.pack(expand=True, fill='both')
 
-        self.entries = Tkinter.Frame(self.topframe)
+        self.entries = tkinter.Frame(self.topframe)
         self.entries.pack(expand=True, fill='both')
 
-        self.figure = Tkinter.Frame(self.topframe)
+        self.figure = tkinter.Frame(self.topframe)
         self.figure.pack(expand=True, fill='both')
 
-        header = Tkinter.Label(self.entries,
+        header = tkinter.Label(self.entries,
             text='Define the transfer function:')
         header.grid(row=0, column=0, padx=20, pady=7)
 
@@ -173,7 +178,7 @@ def __init__(self, parent):
         self.tfi = TFInput(self.entries)
         self.sys = self.tfi.get_tf()
 
-        Tkinter.Button(self.entries, text='Apply', command=self.apply,
+        tkinter.Button(self.entries, text='Apply', command=self.apply,
                width=9).grid(row=0, column=1, rowspan=3, padx=10, pady=5)
 
         self.f_bode = plt.figure(figsize=(4, 4))
@@ -183,25 +188,25 @@ def __init__(self, parent):
 
         self.canvas_pzmap = FigureCanvasTkAgg(self.f_pzmap,
                                                 master=self.figure)
-        self.canvas_pzmap.show()
+        self.canvas_pzmap.draw()
         self.canvas_pzmap.get_tk_widget().grid(row=0, column=0,
                                                 padx=0, pady=0)
 
         self.canvas_bode = FigureCanvasTkAgg(self.f_bode,
                                                 master=self.figure)
-        self.canvas_bode.show()
+        self.canvas_bode.draw()
         self.canvas_bode.get_tk_widget().grid(row=0, column=1,
                                                 padx=0, pady=0)
 
         self.canvas_step = FigureCanvasTkAgg(self.f_step,
                                                 master=self.figure)
-        self.canvas_step.show()
+        self.canvas_step.draw()
         self.canvas_step.get_tk_widget().grid(row=1, column=0,
                                                 padx=0, pady=0)
 
         self.canvas_nyquist = FigureCanvasTkAgg(self.f_nyquist, 
                                                 master=self.figure)
-        self.canvas_nyquist.show()
+        self.canvas_nyquist.draw()
         self.canvas_nyquist.get_tk_widget().grid(row=1, column=1,
                                                     padx=0, pady=0)
 
@@ -300,7 +305,7 @@ def mouse_move(self, event):
                 tfcn = self.tfi.get_tf()
             if (tfcn != None):
                 self.draw_pz(tfcn)
-                self.canvas_pzmap.show()
+                self.canvas_pzmap.draw()
 
     def apply(self):
         """Evaluates the transfer function and produces different plots for
@@ -348,13 +353,13 @@ def redraw(self):
             tvec, yvec = control.matlab.step(self.sys)
             plt.plot(tvec.T, yvec)
         except:
-            print "Error plotting step response"
+            print("Error plotting step response")
         plt.suptitle('Step Response')
 
-        self.canvas_pzmap.show()
-        self.canvas_bode.show()
-        self.canvas_step.show()
-        self.canvas_nyquist.show()
+        self.canvas_pzmap.draw()
+        self.canvas_bode.draw()
+        self.canvas_step.draw()
+        self.canvas_nyquist.draw()
 
 def create_analysis():
     """ Create main object """
@@ -364,7 +369,7 @@ def handler():
         sys.exit()
 
     # Launch a GUI for the Analysis module
-    root = Tkinter.Tk()
+    root = tkinter.Tk()
     root.protocol("WM_DELETE_WINDOW", handler)
     Pmw.initialise(root)    
     root.title('Analysis of Linear Systems')
@@ -372,4 +377,6 @@ def handler():
     root.mainloop()
 
 if __name__ == '__main__':
-    create_analysis()
+    import os
+    if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
+        create_analysis()