diff --git a/.travis.yml b/.travis.yml
index c2c672bf..b9826bab 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -5,7 +5,6 @@ python:
 sudo: false
 
 env:
-    - TEST_DIR=tests/em/examples
     - TEST_DIR=tests/em/fdem/forward
     - TEST_DIR=tests/em/fdem/inverse/derivs
     - TEST_DIR=tests/em/fdem/inverse/adjoint
diff --git a/SimPEG/EM/Examples/__init__.py b/SimPEG/EM/Examples/__init__.py
deleted file mode 100644
index eb36678d..00000000
--- a/SimPEG/EM/Examples/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-import CylInversion
diff --git a/SimPEG/EM/Examples/CylInversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py
similarity index 96%
rename from SimPEG/EM/Examples/CylInversion.py
rename to SimPEG/Examples/EM_FDEM_1D_Inversion.py
index cfcfcfc1..fdf2750e 100644
--- a/SimPEG/EM/Examples/CylInversion.py
+++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py
@@ -4,6 +4,13 @@ from scipy.constants import mu_0
 import matplotlib.pyplot as plt
 
 def run(plotIt=True):
+    """
+        EM: FDEM: 1D: Inversion
+        =======================
+
+        Here we will create and run a FDEM 1D inversion.
+
+    """
 
     cs, ncx, ncz, npad = 5., 25, 15, 15
     hx = [(cs,ncx), (cs,npad,1.3)]
diff --git a/SimPEG/FLOW/Examples/Celia1990.py b/SimPEG/Examples/FLOW_Richards_1D_Celia1990.py
similarity index 54%
rename from SimPEG/FLOW/Examples/Celia1990.py
rename to SimPEG/Examples/FLOW_Richards_1D_Celia1990.py
index 24ae82a6..7c663e6b 100644
--- a/SimPEG/FLOW/Examples/Celia1990.py
+++ b/SimPEG/Examples/FLOW_Richards_1D_Celia1990.py
@@ -3,6 +3,39 @@ from SimPEG.FLOW import Richards
 import matplotlib.pyplot as plt
 
 def run(plotIt=True):
+    """
+        FLOW: Richards: 1D: Celia1990
+        =============================
+
+        There are two different forms of Richards equation that differ
+        on how they deal with the non-linearity in the time-stepping term.
+
+        The most fundamental form, referred to as the
+        'mixed'-form of Richards Equation Celia1990_
+
+        .. math::
+
+            \\frac{\partial \\theta(\psi)}{\partial t} - \\nabla \cdot k(\psi) \\nabla \psi - \\frac{\partial k(\psi)}{\partial z} = 0
+            \quad \psi \in \Omega
+
+        where \\\\(\\\\theta\\\\) is water content, and \\\\(\\\\psi\\\\) is pressure head.
+        This formulation of Richards equation is called the
+        'mixed'-form because the equation is parameterized in \\\\(\\\\psi\\\\)
+        but the time-stepping is in terms of \\\\(\\\\theta\\\\).
+
+        As noted in Celia1990_ the 'head'-based form of Richards
+        equation can be written in the continuous form as:
+
+        .. math::
+
+            \\frac{\partial \\theta}{\partial \psi}\\frac{\partial \psi}{\partial t} - \\nabla \cdot k(\psi) \\nabla \psi - \\frac{\partial k(\psi)}{\partial z} = 0 \quad \psi \in \Omega
+
+        However, it can be shown that this does not conserve mass in the discrete formulation.
+
+        Here we reproduce the results from Celia1990_ demonstrating the head-based formulation and the mixed-formulation.
+
+        .. _Celia1990: http://www.webpages.uidaho.edu/ch/papers/Celia.pdf
+    """
     M = Mesh.TensorMesh([np.ones(40)])
     M.setCellGradBC('dirichlet')
     params = Richards.Empirical.HaverkampParams().celia1990
@@ -47,6 +80,7 @@ def run(plotIt=True):
     plt.xlabel('Depth, cm')
     plt.ylabel('Pressure Head, cm')
     plt.legend(('$\Delta t$ = 10 sec','$\Delta t$ = 30 sec','$\Delta t$ = 120 sec'))
+    plt.show()
 
 if __name__ == '__main__':
     run()
diff --git a/SimPEG/Examples/Forward_BasicDirectCurrent.py b/SimPEG/Examples/Forward_BasicDirectCurrent.py
index a3c9ff97..06c9e79e 100644
--- a/SimPEG/Examples/Forward_BasicDirectCurrent.py
+++ b/SimPEG/Examples/Forward_BasicDirectCurrent.py
@@ -73,5 +73,4 @@ def run(plotIt=True):
 
 
 if __name__ == '__main__':
-    Utils._makeExample(__file__)
     run()
diff --git a/SimPEG/Examples/Inversion_Linear.py b/SimPEG/Examples/Inversion_Linear.py
index 2446bd58..a8a0eddc 100644
--- a/SimPEG/Examples/Inversion_Linear.py
+++ b/SimPEG/Examples/Inversion_Linear.py
@@ -1,29 +1,37 @@
 from SimPEG import *
 
-class LinearSurvey(Survey.BaseSurvey):
-    def projectFields(self, u):
-        return u
-
-class LinearProblem(Problem.BaseProblem):
-    """docstring for LinearProblem"""
-
-    surveyPair = LinearSurvey
-
-    def __init__(self, mesh, G, **kwargs):
-        Problem.BaseProblem.__init__(self, mesh, **kwargs)
-        self.G = G
-
-    def fields(self, m, u=None):
-        return self.G.dot(m)
-
-    def Jvec(self, m, v, u=None):
-        return self.G.dot(v)
-
-    def Jtvec(self, m, v, u=None):
-        return self.G.T.dot(v)
-
 
 def run(N=100, plotIt=True):
+    """
+        Inversion: Linear Problem
+        =========================
+
+        Here we go over the basics of creating a linear problem and inversion.
+
+    """
+
+    class LinearSurvey(Survey.BaseSurvey):
+        def projectFields(self, u):
+            return u
+
+    class LinearProblem(Problem.BaseProblem):
+
+        surveyPair = LinearSurvey
+
+        def __init__(self, mesh, G, **kwargs):
+            Problem.BaseProblem.__init__(self, mesh, **kwargs)
+            self.G = G
+
+        def fields(self, m, u=None):
+            return self.G.dot(m)
+
+        def Jvec(self, m, v, u=None):
+            return self.G.dot(v)
+
+        def Jtvec(self, m, v, u=None):
+            return self.G.T.dot(v)
+
+
     np.random.seed(1)
 
     mesh = Mesh.TensorMesh([N])
@@ -79,5 +87,4 @@ def run(N=100, plotIt=True):
     return prob, survey, mesh, mrec
 
 if __name__ == '__main__':
-    Utils._makeExample(__file__)
     run()
diff --git a/SimPEG/Examples/Mesh_Basic_PlotImage.py b/SimPEG/Examples/Mesh_Basic_PlotImage.py
new file mode 100644
index 00000000..3154f281
--- /dev/null
+++ b/SimPEG/Examples/Mesh_Basic_PlotImage.py
@@ -0,0 +1,46 @@
+from SimPEG import *
+
+def run(plotIt=True):
+    """
+        Mesh: Basic: PlotImage
+        ======================
+
+        You can use M.PlotImage to plot images on all of the Meshes.
+
+
+    """
+    M = Mesh.TensorMesh([32,32])
+    v = Utils.ModelBuilder.randomModel(M.vnC, seed=789)
+    v = Utils.mkvc(v)
+
+    O = Mesh.TreeMesh([32,32])
+    O.refine(1)
+    def function(cell):
+        if (cell.center[0] < 0.75 and cell.center[0] > 0.25 and
+            cell.center[1] < 0.75 and cell.center[1] > 0.25):return 5
+        if (cell.center[0] < 0.9 and cell.center[0] > 0.1 and
+            cell.center[1] < 0.9 and cell.center[1] > 0.1):return 4
+        return 3
+    O.refine(function)
+
+    P = M.getInterpolationMat(O.gridCC, 'CC')
+
+    ov = P * v
+
+    if plotIt:
+        import matplotlib.pyplot as plt
+
+        fig, axes = plt.subplots(1,2,figsize=(10,5))
+
+        out = M.plotImage(v, grid=True, ax=axes[0])
+        cb = plt.colorbar(out[0], ax=axes[0]); cb.set_label("Random Field")
+        axes[0].set_title('TensorMesh')
+
+        out = O.plotImage(ov, grid=True, ax=axes[1], clim=[0,1])
+        cb = plt.colorbar(out[0], ax=axes[1]); cb.set_label("Random Field")
+        axes[1].set_title('TreeMesh')
+
+        plt.show()
+
+if __name__ == '__main__':
+    run()
diff --git a/SimPEG/Examples/Mesh_ThreeMeshes.py b/SimPEG/Examples/Mesh_Basic_Types.py
similarity index 82%
rename from SimPEG/Examples/Mesh_ThreeMeshes.py
rename to SimPEG/Examples/Mesh_Basic_Types.py
index 2fc7fa8d..430fe698 100644
--- a/SimPEG/Examples/Mesh_ThreeMeshes.py
+++ b/SimPEG/Examples/Mesh_Basic_Types.py
@@ -1,6 +1,13 @@
 from SimPEG import *
 
 def run(plotIt=True):
+    """
+        Mesh: Basic: Types
+        ==================
+
+        Here we show SimPEG used to create three different types of meshes.
+
+    """
     sz = [16,16]
     tM = Mesh.TensorMesh(sz)
     qM = Mesh.TreeMesh(sz)
@@ -20,5 +27,4 @@ def run(plotIt=True):
         plt.show()
 
 if __name__ == '__main__':
-    Utils._makeExample(__file__)
     run()
diff --git a/SimPEG/Examples/Mesh_QuadTree_Create.py b/SimPEG/Examples/Mesh_QuadTree_Creation.py
similarity index 50%
rename from SimPEG/Examples/Mesh_QuadTree_Create.py
rename to SimPEG/Examples/Mesh_QuadTree_Creation.py
index f4417855..ede69a63 100644
--- a/SimPEG/Examples/Mesh_QuadTree_Create.py
+++ b/SimPEG/Examples/Mesh_QuadTree_Creation.py
@@ -1,7 +1,18 @@
 from SimPEG import *
 
 def run(plotIt=True):
-    from SimPEG import Mesh, np
+    """
+        Mesh: QuadTree: Creation
+        ========================
+
+        You can give the refine method a function, which is evaluated on every cell
+        of the TreeMesh.
+
+        Occasionally it is useful to initially refine to a constant level
+        (e.g. 3 in this 32x32 mesh). This means the function is first evaluated
+        on an 8x8 mesh (2^3).
+
+    """
     M = Mesh.TreeMesh([32,32])
     M.refine(3)
     def function(cell):
@@ -14,5 +25,4 @@ def run(plotIt=True):
     if plotIt: M.plotGrid(showIt=True)
 
 if __name__ == '__main__':
-    Utils._makeExample(__file__)
     run()
diff --git a/SimPEG/Examples/Mesh_QuadTree_HangingNodes.py b/SimPEG/Examples/Mesh_QuadTree_HangingNodes.py
new file mode 100644
index 00000000..11726995
--- /dev/null
+++ b/SimPEG/Examples/Mesh_QuadTree_HangingNodes.py
@@ -0,0 +1,32 @@
+from SimPEG import *
+
+def run(plotIt=True):
+    """
+        Mesh: QuadTree: Hanging Nodes
+        =============================
+
+        You can give the refine method a function, which is evaluated on every cell
+        of the TreeMesh.
+
+        Occasionally it is useful to initially refine to a constant level
+        (e.g. 3 in this 32x32 mesh). This means the function is first evaluated
+        on an 8x8 mesh (2^3).
+
+    """
+    M = Mesh.TreeMesh([8,8])
+    def function(cell):
+        xyz = cell.center
+        dist = ((xyz - [0.25,0.25])**2).sum()**0.5
+        if dist < 0.25:
+            return 3
+        return 2
+    M.refine(function);
+    M.number()
+    if plotIt:
+        import matplotlib.pyplot as plt
+        M.plotGrid(nodes=True, cells=True, facesX=True)
+        plt.legend(('Grid', 'Cell Centers', 'Nodes', 'Hanging Nodes', 'X faces', 'Hanging X faces'))
+        plt.show()
+
+if __name__ == '__main__':
+    run()
diff --git a/SimPEG/Examples/Mesh_Tensor_Creation.py b/SimPEG/Examples/Mesh_Tensor_Creation.py
new file mode 100644
index 00000000..31ad3d69
--- /dev/null
+++ b/SimPEG/Examples/Mesh_Tensor_Creation.py
@@ -0,0 +1,35 @@
+from SimPEG import *
+
+def run(plotIt=True):
+    """
+
+        Mesh: Tensor: Creation
+        ======================
+
+        For tensor meshes, there are some functions that can come
+        in handy. For example, creating mesh tensors can be a bit time
+        consuming, these can be created speedily by just giving numbers
+        and sizes of padding. See the example below, that follows this
+        notation::
+
+            h1 = (
+                   (cellSize, numPad, [, increaseFactor]),
+                   (cellSize, numCore),
+                   (cellSize, numPad, [, increaseFactor])
+                 )
+
+        .. note::
+
+            You can center your mesh by passing a 'C' for the x0[i] position.
+            A 'N' will make the entire mesh negative, and a '0' (or a 0) will
+            make the mesh start at zero.
+
+    """
+    h1 = [(10, 5, -1.3), (5, 20), (10, 3, 1.3)]
+    M = Mesh.TensorMesh([h1, h1], x0='CN')
+    if plotIt:
+        M.plotGrid(showIt=True)
+
+if __name__ == '__main__':
+    run()
+
diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py
index 53b1f9a8..ac219a7f 100644
--- a/SimPEG/Examples/__init__.py
+++ b/SimPEG/Examples/__init__.py
@@ -5,5 +5,63 @@ __all__ = []
 for x in glob(p.join(p.dirname(__file__), '*.py')):
     if not p.basename(x).startswith('__'):
         __import__(p.basename(x)[:-3], globals(), locals())
-        __all__ += [p.basename(x)]
+        __all__ += [p.basename(x)[:-3]]
 del glob, p, x
+
+if __name__ == '__main__':
+    """
+
+        Run the following to create the examples documentation.
+
+    """
+
+    import shutil, os
+    from SimPEG import Examples
+
+    def _makeExample(filePath, runFunction):
+        filePath = os.path.realpath(filePath)
+        name = filePath.split(os.path.sep)[-1].rstrip('.pyc').rstrip('.py')
+
+        docstr = runFunction.__doc__
+        if docstr is None:
+            doc = '%s\n%s'%(name.replace('_',' '),'='*len(name))
+        else:
+            doc = '\n'.join([_[8:].rstrip() for _ in docstr.split('\n')])
+
+        out = """.. _examples_%s:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+%s
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.%s.run()
+
+.. literalinclude:: ../../SimPEG/Examples/%s.py
+    :language: python
+    :linenos:
+"""%(name,doc,name,name)
+
+        rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst']))
+
+        f = open(rst, 'w')
+        f.write(out)
+        f.close()
+
+
+    docExamplesDir = os.path.sep.join(os.path.realpath(__file__).split(os.path.sep)[:-3] + ['docs', 'examples'])
+    shutil.rmtree(docExamplesDir)
+    os.makedirs(docExamplesDir)
+
+    for ex in dir(Examples):
+        if ex.startswith('_'): continue
+        E = getattr(Examples,ex)
+        _makeExample(E.__file__, E.run)
diff --git a/SimPEG/FLOW/Examples/__init__.py b/SimPEG/FLOW/Examples/__init__.py
deleted file mode 100644
index 7a894e11..00000000
--- a/SimPEG/FLOW/Examples/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-import Celia1990
diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py
index 69a44fdc..263909e4 100644
--- a/SimPEG/Mesh/TreeMesh.py
+++ b/SimPEG/Mesh/TreeMesh.py
@@ -1975,24 +1975,28 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
             fig = ax.figure
 
         if grid:
+            X, Y, Z = [], [], []
             for ind in self._sortedCells:
                 p = self._asPointer(ind)
                 n = self._cellN(p)
                 h = self._cellH(p)
-                x = [n[0]    , n[0] + h[0], n[0] + h[0], n[0]       , n[0]]
-                y = [n[1]    , n[1]       , n[1] + h[1], n[1] + h[1], n[1]]
                 if self.dim == 2:
-                    ax.plot(x,y, 'b-')
+                    X += [n[0]    , n[0] + h[0], n[0] + h[0], n[0]       , n[0], np.nan]
+                    Y += [n[1]    , n[1]       , n[1] + h[1], n[1] + h[1], n[1], np.nan]
                 elif self.dim == 3:
-                    ax.plot(x,y, 'b-', zs=[n[2]]*5)
-                    z = [n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2]]
-                    ax.plot(x,y, 'b-', zs=z)
+                    X += [n[0]    , n[0] + h[0], n[0] + h[0], n[0]       , n[0], np.nan]*2
+                    Y += [n[1]    , n[1]       , n[1] + h[1], n[1] + h[1], n[1], np.nan]*2
+                    Z += [n[2]]*5+[np.nan]
+                    Z += [n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], np.nan]
                     sides = [0,0], [h[0],0], [0,h[1]], [h[0],h[1]]
                     for s in sides:
-                        x = [n[0] + s[0], n[0] + s[0]]
-                        y = [n[1] + s[1], n[1] + s[1]]
-                        z = [n[2]       , n[2] + h[2]]
-                        ax.plot(x,y, 'b-', zs=z)
+                        X += [n[0] + s[0], n[0] + s[0]]
+                        Y += [n[1] + s[1], n[1] + s[1]]
+                        Z += [n[2]       , n[2] + h[2]]
+            if self.dim == 2:
+                ax.plot(X,Y, 'b-')
+            elif self.dim == 3:
+                ax.plot(X,Y, 'b-', zs=Z)
 
         if self.dim == 2:
             if cells:
@@ -2004,11 +2008,11 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
                 ax.plot(self._gridN[:,0], self._gridN[:,1], 'ms')
                 ax.plot(self._gridN[self._hangingN.keys(),0], self._gridN[self._hangingN.keys(),1], 'ms', ms=10, mfc='none', mec='m')
             if facesX:
-                ax.plot(self._gridFx[self._hangingFx.keys(),0], self._gridFx[self._hangingFx.keys(),1], 'gs', ms=10, mfc='none', mec='g')
                 ax.plot(self._gridFx[:,0], self._gridFx[:,1], 'g>')
+                ax.plot(self._gridFx[self._hangingFx.keys(),0], self._gridFx[self._hangingFx.keys(),1], 'gs', ms=10, mfc='none', mec='g')
             if facesY:
-                ax.plot(self._gridFy[self._hangingFy.keys(),0], self._gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g')
                 ax.plot(self._gridFy[:,0], self._gridFy[:,1], 'g^')
+                ax.plot(self._gridFy[self._hangingFy.keys(),0], self._gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g')
             ax.set_xlabel('x1')
             ax.set_ylabel('x2')
         elif self.dim == 3:
@@ -2079,12 +2083,15 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
         ax.grid(True)
         if showIt:plt.show()
 
-    def plotImage(self, I, ax=None, showIt=True, grid=False):
+    def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None):
         if self.dim == 3: raise Exception('Use plot slice?')
 
         if ax is None: ax = plt.subplot(111)
         jet = cm = plt.get_cmap('jet')
-        cNorm  = colors.Normalize(vmin=I.min(), vmax=I.max())
+        cNorm  = colors.Normalize(
+            vmin=I.min() if clim is None else clim[0],
+            vmax=I.max() if clim is None else clim[1])
+
         scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)
         ax.set_xlim((self.x0[0], self.h[0].sum()))
         ax.set_ylim((self.x0[1], self.h[1].sum()))
@@ -2093,8 +2100,10 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
             ax.add_patch(plt.Rectangle((x0[0], x0[1]), sz[0], sz[1], facecolor=scalarMap.to_rgba(I[ii]), edgecolor='k' if grid else 'none'))
             # if text: ax.text(self.center[0],self.center[1],self.num)
         scalarMap._A = []  # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots
-        plt.colorbar(scalarMap)
+        ax.set_xlabel('x')
+        ax.set_ylabel('y')
         if showIt: plt.show()
+        return [scalarMap]
 
     def plotSlice(self, v, vType='CC',
         normal='Z', ind=None, grid=True, view='real',
diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py
index 67e1b117..5280ae79 100644
--- a/SimPEG/Utils/__init__.py
+++ b/SimPEG/Utils/__init__.py
@@ -1,6 +1,5 @@
 from matutils import *
 from codeutils import *
-from codeutils import _makeExample
 from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile
 from curvutils import volTetra, faceInfo, indexCube
 from interputils import interpmat
diff --git a/SimPEG/Utils/codeutils.py b/SimPEG/Utils/codeutils.py
index fef34282..4a9a28a7 100644
--- a/SimPEG/Utils/codeutils.py
+++ b/SimPEG/Utils/codeutils.py
@@ -227,35 +227,3 @@ def requires(var):
 
         return requiresVarWrapper
     return requiresVar
-
-def _makeExample(filePath):
-
-    import os
-    name = filePath.split(os.path.sep)[-1][:-3]
-    out = """.. _examples_%s:
-
-.. ------------------------------ ..
-..  THIS FILE IS AUTO GENEREATED  ..
-.. ------------------------------ ..
-
-%s
-%s
-
-.. plot::
-
-    from SimPEG import Examples
-    Examples.%s.run()
-
-.. literalinclude:: ../../SimPEG/Examples/%s.py
-    :language: python
-    :linenos:
-"""%(name,name.replace('_',' '),'='*len(name),name,name)
-
-    rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst']))
-
-    f = open(rst, 'w')
-    f.write(out)
-    f.close()
-
-
-
diff --git a/docs/examples/EM_FDEM_1D_Inversion.rst b/docs/examples/EM_FDEM_1D_Inversion.rst
new file mode 100644
index 00000000..acbc8cdc
--- /dev/null
+++ b/docs/examples/EM_FDEM_1D_Inversion.rst
@@ -0,0 +1,26 @@
+.. _examples_EM_FDEM_1D_Inversion:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+EM: FDEM: 1D: Inversion
+=======================
+
+Here we will create and run a FDEM 1D inversion.
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.EM_FDEM_1D_Inversion.run()
+
+.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_1D_Inversion.py
+    :language: python
+    :linenos:
diff --git a/docs/examples/FLOW_Richards_1D_Celia1990.rst b/docs/examples/FLOW_Richards_1D_Celia1990.rst
new file mode 100644
index 00000000..d2e01c13
--- /dev/null
+++ b/docs/examples/FLOW_Richards_1D_Celia1990.rst
@@ -0,0 +1,52 @@
+.. _examples_FLOW_Richards_1D_Celia1990:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+FLOW: Richards: 1D: Celia1990
+=============================
+
+There are two different forms of Richards equation that differ
+on how they deal with the non-linearity in the time-stepping term.
+
+The most fundamental form, referred to as the
+'mixed'-form of Richards Equation Celia1990_
+
+.. math::
+
+    \frac{\partial \theta(\psi)}{\partial t} - \nabla \cdot k(\psi) \nabla \psi - \frac{\partial k(\psi)}{\partial z} = 0
+    \quad \psi \in \Omega
+
+where \\(\\theta\\) is water content, and \\(\\psi\\) is pressure head.
+This formulation of Richards equation is called the
+'mixed'-form because the equation is parameterized in \\(\\psi\\)
+but the time-stepping is in terms of \\(\\theta\\).
+
+As noted in Celia1990_ the 'head'-based form of Richards
+equation can be written in the continuous form as:
+
+.. math::
+
+    \frac{\partial \theta}{\partial \psi}\frac{\partial \psi}{\partial t} - \nabla \cdot k(\psi) \nabla \psi - \frac{\partial k(\psi)}{\partial z} = 0 \quad \psi \in \Omega
+
+However, it can be shown that this does not conserve mass in the discrete formulation.
+
+Here we reproduce the results from Celia1990_ demonstrating the head-based formulation and the mixed-formulation.
+
+.. _Celia1990: http://www.webpages.uidaho.edu/ch/papers/Celia.pdf
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.FLOW_Richards_1D_Celia1990.run()
+
+.. literalinclude:: ../../SimPEG/Examples/FLOW_Richards_1D_Celia1990.py
+    :language: python
+    :linenos:
diff --git a/docs/examples/Forward_BasicDirectCurrent.rst b/docs/examples/Forward_BasicDirectCurrent.rst
index 6fc816e7..20b39eb8 100644
--- a/docs/examples/Forward_BasicDirectCurrent.rst
+++ b/docs/examples/Forward_BasicDirectCurrent.rst
@@ -1,8 +1,12 @@
 .. _examples_Forward_BasicDirectCurrent:
 
-.. ------------------------------ ..
-..  THIS FILE IS AUTO GENEREATED  ..
-.. ------------------------------ ..
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
 
 Forward BasicDirectCurrent
 ==========================
diff --git a/docs/examples/Inversion_Linear.rst b/docs/examples/Inversion_Linear.rst
index 637315ae..d635d8e1 100644
--- a/docs/examples/Inversion_Linear.rst
+++ b/docs/examples/Inversion_Linear.rst
@@ -1,11 +1,20 @@
 .. _examples_Inversion_Linear:
 
-.. ------------------------------ ..
-..  THIS FILE IS AUTO GENEREATED  ..
-.. ------------------------------ ..
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+Inversion: Linear Problem
+=========================
+
+Here we go over the basics of creating a linear problem and inversion.
+
 
-Inversion Linear
-================
 
 .. plot::
 
diff --git a/docs/examples/Mesh_Basic_PlotImage.rst b/docs/examples/Mesh_Basic_PlotImage.rst
new file mode 100644
index 00000000..a730f303
--- /dev/null
+++ b/docs/examples/Mesh_Basic_PlotImage.rst
@@ -0,0 +1,27 @@
+.. _examples_Mesh_Basic_PlotImage:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+Mesh: Basic: PlotImage
+======================
+
+You can use M.PlotImage to plot images on all of the Meshes.
+
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.Mesh_Basic_PlotImage.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Mesh_Basic_PlotImage.py
+    :language: python
+    :linenos:
diff --git a/docs/examples/Mesh_Basic_Types.rst b/docs/examples/Mesh_Basic_Types.rst
new file mode 100644
index 00000000..9bbce0e8
--- /dev/null
+++ b/docs/examples/Mesh_Basic_Types.rst
@@ -0,0 +1,26 @@
+.. _examples_Mesh_Basic_Types:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+Mesh: Basic: Types
+==================
+
+Here we show SimPEG used to create three different types of meshes.
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.Mesh_Basic_Types.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Mesh_Basic_Types.py
+    :language: python
+    :linenos:
diff --git a/docs/examples/Mesh_QuadTree_Create.rst b/docs/examples/Mesh_QuadTree_Create.rst
deleted file mode 100644
index 644593d5..00000000
--- a/docs/examples/Mesh_QuadTree_Create.rst
+++ /dev/null
@@ -1,17 +0,0 @@
-.. _examples_Mesh_QuadTree_Create:
-
-.. ------------------------------ ..
-..  THIS FILE IS AUTO GENEREATED  ..
-.. ------------------------------ ..
-
-Mesh QuadTree Create
-====================
-
-.. plot::
-
-    from SimPEG import Examples
-    Examples.Mesh_QuadTree_Create.run()
-
-.. literalinclude:: ../../SimPEG/Examples/Mesh_QuadTree_Create.py
-    :language: python
-    :linenos:
diff --git a/docs/examples/Mesh_QuadTree_Creation.rst b/docs/examples/Mesh_QuadTree_Creation.rst
new file mode 100644
index 00000000..5db5a982
--- /dev/null
+++ b/docs/examples/Mesh_QuadTree_Creation.rst
@@ -0,0 +1,31 @@
+.. _examples_Mesh_QuadTree_Creation:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+Mesh: QuadTree: Creation
+========================
+
+You can give the refine method a function, which is evaluated on every cell
+of the TreeMesh.
+
+Occasionally it is useful to initially refine to a constant level
+(e.g. 3 in this 32x32 mesh). This means the function is first evaluated
+on an 8x8 mesh (2^3).
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.Mesh_QuadTree_Creation.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Mesh_QuadTree_Creation.py
+    :language: python
+    :linenos:
diff --git a/docs/examples/Mesh_QuadTree_HangingNodes.rst b/docs/examples/Mesh_QuadTree_HangingNodes.rst
new file mode 100644
index 00000000..93875478
--- /dev/null
+++ b/docs/examples/Mesh_QuadTree_HangingNodes.rst
@@ -0,0 +1,31 @@
+.. _examples_Mesh_QuadTree_HangingNodes:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+Mesh: QuadTree: Hanging Nodes
+=============================
+
+You can give the refine method a function, which is evaluated on every cell
+of the TreeMesh.
+
+Occasionally it is useful to initially refine to a constant level
+(e.g. 3 in this 32x32 mesh). This means the function is first evaluated
+on an 8x8 mesh (2^3).
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.Mesh_QuadTree_HangingNodes.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Mesh_QuadTree_HangingNodes.py
+    :language: python
+    :linenos:
diff --git a/docs/examples/Mesh_Tensor_Creation.rst b/docs/examples/Mesh_Tensor_Creation.rst
new file mode 100644
index 00000000..e6cc5b67
--- /dev/null
+++ b/docs/examples/Mesh_Tensor_Creation.rst
@@ -0,0 +1,43 @@
+.. _examples_Mesh_Tensor_Creation:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+
+Mesh: Tensor: Creation
+======================
+
+For tensor meshes, there are some functions that can come
+in handy. For example, creating mesh tensors can be a bit time
+consuming, these can be created speedily by just giving numbers
+and sizes of padding. See the example below, that follows this
+notation::
+
+    h1 = (
+           (cellSize, numPad, [, increaseFactor]),
+           (cellSize, numCore),
+           (cellSize, numPad, [, increaseFactor])
+         )
+
+.. note::
+
+    You can center your mesh by passing a 'C' for the x0[i] position.
+    A 'N' will make the entire mesh negative, and a '0' (or a 0) will
+    make the mesh start at zero.
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.Mesh_Tensor_Creation.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Mesh_Tensor_Creation.py
+    :language: python
+    :linenos:
diff --git a/docs/examples/Mesh_ThreeMeshes.rst b/docs/examples/Mesh_ThreeMeshes.rst
deleted file mode 100644
index 948398a7..00000000
--- a/docs/examples/Mesh_ThreeMeshes.rst
+++ /dev/null
@@ -1,17 +0,0 @@
-.. _examples_Mesh_ThreeMeshes:
-
-.. ------------------------------ ..
-..  THIS FILE IS AUTO GENEREATED  ..
-.. ------------------------------ ..
-
-Mesh ThreeMeshes
-================
-
-.. plot::
-
-    from SimPEG import Examples
-    Examples.Mesh_ThreeMeshes.run()
-
-.. literalinclude:: ../../SimPEG/Examples/Mesh_ThreeMeshes.py
-    :language: python
-    :linenos:
diff --git a/tests/em/examples/__init__.py b/tests/em/examples/__init__.py
deleted file mode 100644
index 38d84328..00000000
--- a/tests/em/examples/__init__.py
+++ /dev/null
@@ -1,11 +0,0 @@
-if __name__ == '__main__':
-    import os
-    import glob
-    import unittest
-    test_file_strings = glob.glob('test_*.py')
-    module_strings = [str[0:len(str)-3] for str in test_file_strings]
-    suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
-              in module_strings]
-    testSuite = unittest.TestSuite(suites)
-
-    unittest.TextTestRunner(verbosity=2).run(testSuite)
diff --git a/tests/em/examples/test_Examples.py b/tests/em/examples/test_Examples.py
deleted file mode 100644
index 5a601d3b..00000000
--- a/tests/em/examples/test_Examples.py
+++ /dev/null
@@ -1,10 +0,0 @@
-import unittest, os
-from SimPEG.EM import Examples
-
-class EM_ExamplesRunning(unittest.TestCase):
-
-    def test_CylInversion(self):
-        Examples.CylInversion.run(plotIt=False)
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/tests/flow/test_examples.py b/tests/flow/test_examples.py
deleted file mode 100644
index bc519e6d..00000000
--- a/tests/flow/test_examples.py
+++ /dev/null
@@ -1,12 +0,0 @@
-import unittest
-import sys
-from SimPEG.FLOW.Examples import Celia1990
-import numpy as np
-
-class TestCelia1990(unittest.TestCase):
-    def test_running(self):
-        Celia1990.run(plotIt=False)
-        self.assertTrue(True)
-
-if __name__ == '__main__':
-    unittest.main()