diff --git a/wafo/covariance/core.py b/wafo/covariance/core.py
index afb88b2..7a8c67e 100644
--- a/wafo/covariance/core.py
+++ b/wafo/covariance/core.py
@@ -68,7 +68,7 @@ def rndnormnd(mean, cov, cases=1):
     Example
     -------
     >>> mu = [0, 5]
-    >>> S = [[1 0.45], [0.45 0.25]]
+    >>> S = [[1 0.45], [0.45, 0.25]]
     >>> r = rndnormnd(mu, S, 1)
 
     plot(r(:,1),r(:,2),'.')
@@ -208,7 +208,8 @@ class CovData1D(PlotData):
         >>> S = Sj.tospecdata()
         >>> R2 = S.tocovdata()
         >>> S1 = R2.tospecdata()
-        >>> abs(S1.data-S.data).max()
+        >>> abs(S1.data-S.data).max() < 1e-4
+        True
 
         >>> S1.plot('r-')
         >>> S.plot('b:')
diff --git a/wafo/covariance/estimation.py b/wafo/covariance/estimation.py
index 51a2ea2..1a2a195 100644
--- a/wafo/covariance/estimation.py
+++ b/wafo/covariance/estimation.py
@@ -110,7 +110,8 @@ class CovarianceEstimator(object):
          >>> x = wafo.data.sea()
          >>> ts = wo.mat2timeseries(x)
          >>> acf = ts.tocovdata(150)
-         >>> h = acf.plot()
+
+         h = acf.plot()
         '''
         lag = self.lag
         window = self.window
diff --git a/wafo/data/info.py b/wafo/data/info.py
index 68a4e79..6eb1206 100644
--- a/wafo/data/info.py
+++ b/wafo/data/info.py
@@ -20,7 +20,7 @@ included in WAFO
 
 """
 import numpy as np
-from numpy import (loadtxt, nan)
+
 import os
 __path2data = os.path.dirname(os.path.realpath(__file__))
 
@@ -31,7 +31,7 @@ _NANS = set(['nan', 'NaN', '-1.#IND00+00', '1.#IND00+00', '-1.#INF00+00'])
 
 
 def _tofloat(x):
-    return nan if x in _NANS else float(x or 0)
+    return np.nan if x in _NANS else float(x or 0)
 
 
 _MYCONVERTER = {}
@@ -42,13 +42,13 @@ for i in range(2):
 def _load(file):  # @ReservedAssignment
     """ local load function
     """
-    return loadtxt(os.path.join(__path2data, file))
+    return np.loadtxt(os.path.join(__path2data, file))
 
 
 def _loadnan(file):  # @ReservedAssignment
     """ local load function accepting nan's
     """
-    return loadtxt(os.path.join(__path2data, file), converters=_MYCONVERTER)
+    return np.loadtxt(os.path.join(__path2data, file), converters=_MYCONVERTER)
 
 
 def atlantic():
@@ -188,8 +188,8 @@ def gfaksr89():
     >>> import pylab
     >>> import wafo
     >>> x = wafo.data.gfaksr89()
-    >>> h = pylab.plot(x[:,0],x[:,1])
 
+    h = pylab.plot(x[:,0],x[:,1])
 
     See also
     --------
@@ -378,7 +378,8 @@ def sfa89():
     >>> import pylab
     >>> import wafo
     >>> x = wafo.data.sfa89()
-    >>> h = pylab.plot(x[:,0],x[:,1])
+
+    h = pylab.plot(x[:,0],x[:,1])
 
     Acknowledgement
     ----------------
@@ -415,7 +416,8 @@ def sn():
     >>> import pylab
     >>> import wafo
     >>> x = wafo.data.sn()
-    >>> h = pylab.plot(x[:,0],x[:,1])
+
+    h = pylab.plot(x[:,0],x[:,1])
 
     See also
     --------
@@ -463,7 +465,8 @@ def yura87():
     >>> import pylab
     >>> import wafo
     >>> x = wafo.data.yura87()
-    >>> h = pylab.plot(x[:,0],x[:,1])
+
+    h = pylab.plot(x[:,0],x[:,1])
 
     Acknowledgement:
     -----------------
diff --git a/wafo/spectrum/core.py b/wafo/spectrum/core.py
index 97e8656..609ca4b 100644
--- a/wafo/spectrum/core.py
+++ b/wafo/spectrum/core.py
@@ -1072,6 +1072,7 @@ class SpecData1D(PlotData):
         >>> w = np.linspace(0,4,256)
         >>> S1 = Sj.tospecdata(w)   #Make spectrum object from numerical values
         >>> S = sm.SpecData1D(Sj(w),w) # Alternatively do it manually
+
         mm = S.to_mm_pdf()
         mm.plot()
         mm.plot(plotflag=1)
@@ -1174,12 +1175,11 @@ class SpecData1D(PlotData):
         >>> Sj = sm.Jonswap()
         >>> S = Sj.tospecdata()
         >>> f = S.to_t_pdf(pdef='Tc', paramt=(0, 10, 51), speed=7)
-        >>> h = f.plot()
-
-        estimated error bounds
-        >>> h2 = plb.plot(f.args, f.data+f.err, 'r', f.args, f.data-f.err, 'r')
 
-        >>> plb.close('all')
+        h = f.plot()
+        # estimated error bounds
+        h2 = plb.plot(f.args, f.data+f.err, 'r', f.args, f.data-f.err, 'r')
+        plb.close('all')
 
         See also
         --------
@@ -3955,11 +3955,13 @@ class SpecData2D(PlotData):
         >>> D = sm.Spreading()
         >>> SD = D.tospecdata2d(sm.Jonswap().tospecdata(),nt=101)
         >>> m,mtext = SD.moment(nr=2,vari='xyt')
-        >>> np.round(m,3),mtext
-        (array([ 3.061,  0.132, -0.   ,  2.13 ,  0.011,  0.008,  1.677, -0.,
-                0.109,  0.109]),
-                ['m0', 'mx', 'my', 'mt', 'mxx', 'myy', 'mtt', 'mxy', 'mxt',
-                'myt'])
+        >>> np.allclose(np.round(m,3),
+        ... [ 3.061,  0.132, -0.   ,  2.13 ,  0.011,  0.008,  1.677, -0.,
+        ...     0.109,  0.109])
+        True
+        >>> mtext == ['m0', 'mx', 'my', 'mt', 'mxx', 'myy', 'mtt', 'mxy',
+        ...            'mxt', 'myt']
+        True
 
         References
         ----------
diff --git a/wafo/spectrum/models.py b/wafo/spectrum/models.py
index 8a9383e..888fa66 100644
--- a/wafo/spectrum/models.py
+++ b/wafo/spectrum/models.py
@@ -297,17 +297,18 @@ def jonswap_peakfact(Hm0, Tp):
     >>> [T,H] = plb.meshgrid(Tp,Hm0)
     >>> gam = wsm.jonswap_peakfact(H,T)
     >>> v = plb.arange(0,8)
-    >>> h = plb.contourf(Tp,Hm0,gam,v);h=plb.colorbar()
 
     >>> Hm0 = plb.arange(1,11)
     >>> Tp  = plb.linspace(2,16)
     >>> T,H = plb.meshgrid(Tp,Hm0)
     >>> gam = wsm.jonswap_peakfact(H,T)
-    >>> h = plb.plot(Tp,gam.T)
-    >>> h = plb.xlabel('Tp [s]')
-    >>> h = plb.ylabel('Peakedness parameter')
 
-    >>> plb.close('all')
+    h = plb.contourf(Tp,Hm0,gam,v);h=plb.colorbar()
+    h = plb.plot(Tp,gam.T)
+    h = plb.xlabel('Tp [s]')
+    h = plb.ylabel('Peakedness parameter')
+
+    plb.close('all')
 
     See also
     --------
@@ -495,13 +496,13 @@ class Jonswap(ModelSpectrum):
     >>> import pylab as plb
     >>> import wafo.spectrum.models as wsm
     >>> S = wsm.Jonswap(Hm0=7, Tp=11,gamma=1)
-    >>> w = plb.linspace(0,5)
-    >>> h = plb.plot(w,S(w))
-
     >>> S2 = wsm.Bretschneider(Hm0=7, Tp=11)
+    >>> w = plb.linspace(0,5)
     >>> all(abs(S(w)-S2(w))<1.e-7)
     True
-    >>> plb.close('all')
+
+    h = plb.plot(w,S(w))
+    plb.close('all')
 
     See also
     --------
@@ -755,11 +756,12 @@ class Tmaspec(Jonswap):
     >>> import pylab as plb
     >>> w = plb.linspace(0,2.5)
     >>> S = wsm.Tmaspec(h=10,gamma=1) # Bretschneider spectrum Hm0=7, Tp=11
-    >>> o=plb.plot(w,S(w))
-    >>> o=plb.plot(w,S(w,h=21))
-    >>> o=plb.plot(w,S(w,h=42))
-    >>> plb.show()
-    >>> plb.close('all')
+
+    o=plb.plot(w,S(w))
+    o=plb.plot(w,S(w,h=21))
+    o=plb.plot(w,S(w,h=42))
+    plb.show()
+    plb.close('all')
 
     See also
     ---------
@@ -850,7 +852,8 @@ class Torsethaugen(ModelSpectrum):
     >>> import pylab as plb
     >>> w = plb.linspace(0,4)
     >>> S = wsm.Torsethaugen(Hm0=6, Tp=8)
-    >>> h=plb.plot(w,S(w),w,S.wind(w),w,S.swell(w))
+
+    h=plb.plot(w,S(w),w,S.wind(w),w,S.swell(w))
 
     See also
     --------
@@ -1429,25 +1432,26 @@ class Spreading(object):
     # Make directionale spectrum
     >>> S = wsm.Jonswap().tospecdata()
     >>> SD = D.tospecdata2d(S)
-    >>> h = SD.plot()
 
     >>> w = plb.linspace(0,3,257)
     >>> theta = plb.linspace(-pi,pi,129)
-    >>> t = plb.contour(D(theta,w)[0].squeeze())
 
     # Make frequency dependent direction spreading
     >>> theta0 = lambda w: w*plb.pi/6.0
     >>> D2 = wsm.Spreading('cos2s',theta0=theta0)
-    >>> t = plb.contour(D2(theta,w)[0])
+
+    h = SD.plot()
+    t = plb.contour(D(theta,w)[0].squeeze())
+    t = plb.contour(D2(theta,w)[0])
 
     # Plot all spreading functions
-    >>> alltypes = ('cos2s','box','mises','poisson','sech2','wrap_norm')
-    >>> for ix in range(len(alltypes)):
+    alltypes = ('cos2s','box','mises','poisson','sech2','wrap_norm')
+    for ix in range(len(alltypes)):
     ...     D3 = wsm.Spreading(alltypes[ix])
     ...     t = plb.figure(ix)
     ...     t = plb.contour(D3(theta,w)[0])
     ...     t = plb.title(alltypes[ix])
-    >>> plb.close('all')
+    plb.close('all')
 
 
     See also
@@ -2007,7 +2011,8 @@ class Spreading(object):
          >>> S = wsm.Jonswap().tospecdata()
          >>> D = wsm.Spreading('cos2s')
          >>> SD = D.tospecdata2d(S)
-         >>> h = SD.plot()
+
+         h = SD.plot()
 
          See also  spreading, rotspec, jonswap, torsethaugen
         '''
diff --git a/wafo/transform/estimation.py b/wafo/transform/estimation.py
index b282d5c..8d85e93 100644
--- a/wafo/transform/estimation.py
+++ b/wafo/transform/estimation.py
@@ -196,7 +196,7 @@ class TransformEstimator(object):
         >>> tp = ts.turning_points()
         >>> mm = tp.cycle_pairs()
         >>> lc = mm.level_crossings()
-        >>> g0, g0emp = lc.trdata(monitor=True) # Monitor the development
+        >>> g0, g0emp = lc.trdata(monitor=False) # Monitor the development
         >>> g1, g1emp = lc.trdata(gvar=0.5 ) # Equal weight on all points
         >>> g2, g2emp = lc.trdata(gvar=[3.5, 0.5, 3.5])  # Less weight on ends
         >>> int(S.tr.dist2gauss()*100)
@@ -213,6 +213,7 @@ class TransformEstimator(object):
         g0.plot() # Check the fit.
 
         See also
+        --------
           troptset, dat2tr, trplot, findcross, smooth
 
         NB! the transformated data will be N(0,1)
@@ -389,7 +390,7 @@ class TransformEstimator(object):
         ...        sigma=Hs/4, ysigma=Hs/4)
         >>> xs = S.sim(ns=2**16, iseed=10)
         >>> ts = mat2timeseries(xs)
-        >>> g0, g0emp = ts.trdata(monitor=True)
+        >>> g0, g0emp = ts.trdata(monitor=False)
         >>> g1, g1emp = ts.trdata(method='m', gvar=0.5 )
         >>> g2, g2emp = ts.trdata(method='n', gvar=[3.5, 0.5, 3.5])
         >>> int(S.tr.dist2gauss()*100)
diff --git a/wafo/wave_theory/core.py b/wafo/wave_theory/core.py
index 297f7c5..2da24c5 100644
--- a/wafo/wave_theory/core.py
+++ b/wafo/wave_theory/core.py
@@ -245,14 +245,15 @@ class TransferFunction(object):
     >>> stypes = ['n', 'n_x', 'n_y'];
     >>> tf = TransferFunction(pos=(0, 0, 0), h=50)
     >>> vals = []
-    >>> fh = plt.plot(t, eta0.real, 'r.')
-    >>> plt.hold(True)
     >>> for i,stype in enumerate(stypes):
     ...    tf.sensortype = stype
     ...    Hw, Gwt = tf.tran(w0,th0)
     ...    vals.append((Hw*Gwt*eta0).real.ravel())
 
-    fh = plt.plot(t, vals[i])
+    fh = plt.plot(t, eta0.real, 'r.')
+    plt.hold(True)
+    for val in vals:
+        fh = plt.plot(t, val)
     plt.show()
 
 
@@ -543,8 +544,10 @@ def wave_pressure(z, Hm0, h=10000, g=9.81, rho=1028):
     -----
     >>> import pylab as plt
     >>> z = -np.linspace(10,20)
-    >>> fh = plt.plot(z, wave_pressure(z, Hm0=1, h=20))
-    >>> plt.show()
+    >>> p = wave_pressure(z, Hm0=1, h=20)
+
+    fh = plt.plot(z, p)
+    plt.show()
 
     See also
     --------
diff --git a/wafo/wave_theory/dispersion_relation.py b/wafo/wave_theory/dispersion_relation.py
index f05be49..28d30d2 100644
--- a/wafo/wave_theory/dispersion_relation.py
+++ b/wafo/wave_theory/dispersion_relation.py
@@ -5,7 +5,7 @@ k2w - Translates from wave number to frequency
 w2k - Translates from frequency to wave number
 """
 import warnings
-#import numpy as np
+# import numpy as np
 from numpy import (atleast_1d, sqrt, ones_like, zeros_like, arctan2, where,
                    tanh, any, sin, cos, sign, inf,
                    flatnonzero, finfo, cosh, abs)
@@ -124,13 +124,13 @@ def w2k(w, theta=0.0, h=inf, g=9.81, count_limit=100):
     >>> import pylab as plb
     >>> import wafo.wave_theory.dispersion_relation as wsd
     >>> w = plb.linspace(0,3);
-    >>> h = plb.plot(w,w2k(w)[0])
     >>> wsd.w2k(range(4))[0]
     array([ 0.        ,  0.1019368 ,  0.4077472 ,  0.91743119])
     >>> wsd.w2k(range(4),h=20)[0]
     array([ 0.        ,  0.10503601,  0.40774726,  0.91743119])
 
-    >>> plb.close('all')
+    h = plb.plot(w,w2k(w)[0])
+    plb.close('all')
 
     See also
     --------