diff --git a/docs/api_DC.rst b/docs/api_DC.rst
index 5aeeb6d7..b854e0c3 100644
--- a/docs/api_DC.rst
+++ b/docs/api_DC.rst
@@ -129,7 +129,11 @@ where
 
 Here \\(\\bm\\) indicates model parameters in discretized space.
 
-.. plot :: ../simpegDC/Examples/Verification.py
+.. plot::
+
+    import simpegDC as DC
+    DC.Examples.Verification.run(plotIt=True)
+
 
 .. automodule:: simpegDC.DC
     :show-inheritance:
diff --git a/simpegDC/BaseDC.py b/simpegDC/BaseDC.py
index 39037f00..b3380741 100644
--- a/simpegDC/BaseDC.py
+++ b/simpegDC/BaseDC.py
@@ -89,7 +89,7 @@ class ProblemDC(Problem.BaseProblem):
     """
 
     surveyPair = SurveyDC
-    Solver = Solver
+    Solver     = Solver
 
     def __init__(self, mesh, **kwargs):
         Problem.BaseProblem.__init__(self, mesh)
@@ -142,7 +142,7 @@ class ProblemDC(Problem.BaseProblem):
     def fields(self, m):
         self.curModel = m
         A    = self.A
-        Ainv = self.Solver(A)
+        Ainv = self.Solver(A, **self.solverOpts)
         Q    = self.survey.getRhs(self.mesh)
         Phi  = Ainv * Q
         return Phi
@@ -218,7 +218,7 @@ class ProblemDC(Problem.BaseProblem):
         mT_dm = self.mapping.deriv(m)
 
         dCdu = A.T
-        Ainv = self.Solver(dCdu)
+        Ainv = self.Solver(dCdu, **self.solverOpts)
         w = Ainv * PT_x_v
 
         Jtv = 0
diff --git a/simpegDC/Examples/Verification.py b/simpegDC/Examples/Verification.py
index 9ac3ffab..5b722da3 100644
--- a/simpegDC/Examples/Verification.py
+++ b/simpegDC/Examples/Verification.py
@@ -1,51 +1,63 @@
 from SimPEG import *
 import simpegDC as DC
 import matplotlib.pyplot as plt
-cs = 25.
-hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
-hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
-hz = [(cs,7, -1.3),(cs,20)]
-mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN')
-sighalf = 1e-2
-sigma = np.ones(mesh.nC)*sighalf
-xtemp = np.linspace(-150, 150, 21)
-ytemp = np.linspace(-150, 150, 21)
-xyz_rxP = Utils.ndgrid(xtemp-10., ytemp, np.r_[0.])
-xyz_rxN = Utils.ndgrid(xtemp+10., ytemp, np.r_[0.])
-xyz_rxM = Utils.ndgrid(xtemp, ytemp, np.r_[0.])
-fig, ax = plt.subplots(1,1, figsize = (5,5))
-mesh.plotSlice(sigma, grid=True, ax = ax)
-ax.plot(xyz_rxP[:,0],xyz_rxP[:,1], 'w.')
-ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)
-rx = DC.DipoleRx(xyz_rxP, xyz_rxN)
-tx = DC.DipoleTx([-200, 0, -12.5],[+200, 0, -12.5], [rx])
-survey = DC.SurveyDC([tx])
-problem = DC.ProblemDC(mesh)
-problem.pair(survey)
-data = survey.dpred(sigma)
 
-def DChalf(txlocP, txlocN, rxloc, sigma, I=1.):
-    rp = (txlocP.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
-    rn = (txlocN.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
-    rP = np.sqrt(((rxloc-rp)**2).sum(axis=1))
-    rN = np.sqrt(((rxloc-rn)**2).sum(axis=1))
-    return I/(sigma*2.*np.pi)*(1/rP-1/rN)
+def run(plotIt=False):
+    cs = 25.
+    hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
+    hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
+    hz = [(cs,7, -1.3),(cs,20)]
+    mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN')
+    sighalf = 1e-2
+    sigma = np.ones(mesh.nC)*sighalf
+    xtemp = np.linspace(-150, 150, 21)
+    ytemp = np.linspace(-150, 150, 21)
+    xyz_rxP = Utils.ndgrid(xtemp-10., ytemp, np.r_[0.])
+    xyz_rxN = Utils.ndgrid(xtemp+10., ytemp, np.r_[0.])
+    xyz_rxM = Utils.ndgrid(xtemp, ytemp, np.r_[0.])
 
-data_analP = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxP, sighalf)
-data_analN = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxN, sighalf)
-data_anal = data_analP-data_analN
-Data_anal = data_anal.reshape((21, 21), order = 'F')
-Data = data.reshape((21, 21), order = 'F')
-X = xyz_rxM[:,0].reshape((21, 21), order = 'F')
-Y = xyz_rxM[:,1].reshape((21, 21), order = 'F')
+    # if plotIt:
+    #     fig, ax = plt.subplots(1,1, figsize = (5,5))
+    #     mesh.plotSlice(sigma, grid=True, ax = ax)
+    #     ax.plot(xyz_rxP[:,0],xyz_rxP[:,1], 'w.')
+    #     ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)
 
-fig, ax = plt.subplots(1,2, figsize = (12, 5))
-vmin = np.r_[data, data_anal].min()
-vmax = np.r_[data, data_anal].max()
-dat1 = ax[1].contourf(X, Y, Data, 60, vmin = vmin, vmax = vmax)
-dat0 = ax[0].contourf(X, Y, Data_anal, 60, vmin = vmin, vmax = vmax)
-cb0 = plt.colorbar(dat1, orientation = 'horizontal', ax = ax[0])
-cb1 = plt.colorbar(dat1, orientation = 'horizontal', ax = ax[1])
-ax[1].set_title('Analytic')
-ax[0].set_title('Computed')
-plt.show()
+    rx = DC.DipoleRx(xyz_rxP, xyz_rxN)
+    tx = DC.DipoleTx([-200, 0, -12.5],[+200, 0, -12.5], [rx])
+    survey = DC.SurveyDC([tx])
+    problem = DC.ProblemDC(mesh)
+    problem.pair(survey)
+    data = survey.dpred(sigma)
+
+    def DChalf(txlocP, txlocN, rxloc, sigma, I=1.):
+        rp = (txlocP.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
+        rn = (txlocN.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
+        rP = np.sqrt(((rxloc-rp)**2).sum(axis=1))
+        rN = np.sqrt(((rxloc-rn)**2).sum(axis=1))
+        return I/(sigma*2.*np.pi)*(1/rP-1/rN)
+
+    data_analP = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxP, sighalf)
+    data_analN = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxN, sighalf)
+    data_anal = data_analP-data_analN
+    Data_anal = data_anal.reshape((21, 21), order = 'F')
+    Data = data.reshape((21, 21), order = 'F')
+    X = xyz_rxM[:,0].reshape((21, 21), order = 'F')
+    Y = xyz_rxM[:,1].reshape((21, 21), order = 'F')
+
+    if plotIt:
+        fig, ax = plt.subplots(1,2, figsize = (12, 5))
+        vmin = np.r_[data, data_anal].min()
+        vmax = np.r_[data, data_anal].max()
+        dat1 = ax[1].contourf(X, Y, Data, 60, vmin = vmin, vmax = vmax)
+        dat0 = ax[0].contourf(X, Y, Data_anal, 60, vmin = vmin, vmax = vmax)
+        cb0 = plt.colorbar(dat1, orientation = 'horizontal', ax = ax[0])
+        cb1 = plt.colorbar(dat1, orientation = 'horizontal', ax = ax[1])
+        ax[1].set_title('Analytic')
+        ax[0].set_title('Computed')
+        plt.show()
+
+    return np.linalg.norm(data-data_anal)/np.linalg.norm(data_anal)
+
+
+if __name__ == '__main__':
+    print run(plotIt=True)
diff --git a/simpegDC/Examples/WennerArray.py b/simpegDC/Examples/WennerArray.py
index ed37ee68..6238f7d3 100644
--- a/simpegDC/Examples/WennerArray.py
+++ b/simpegDC/Examples/WennerArray.py
@@ -6,6 +6,7 @@ import matplotlib.pyplot as plt
 def getTxList(nElecs, aSpacing, in2D=False, plotIt=False):
 
     elocs = np.arange(0,aSpacing*nElecs,aSpacing)
+    elocs -= (nElecs*aSpacing - aSpacing)/2
     space = 1
     WENNER = np.zeros((0,),dtype=int)
     for ii in range(nElecs):
diff --git a/simpegDC/Examples/__init__.py b/simpegDC/Examples/__init__.py
index c1a78483..99ddc1df 100644
--- a/simpegDC/Examples/__init__.py
+++ b/simpegDC/Examples/__init__.py
@@ -1 +1,2 @@
 import WennerArray
+import Verification
diff --git a/simpegDC/Tests/test_forward_DCproblem_analytic.py b/simpegDC/Tests/test_forward_DCproblem_analytic.py
new file mode 100644
index 00000000..d4677bef
--- /dev/null
+++ b/simpegDC/Tests/test_forward_DCproblem_analytic.py
@@ -0,0 +1,13 @@
+import unittest
+import simpegDC as DC
+
+
+class DCAnalyticTests(unittest.TestCase):
+
+    def test_forwardAnalytic(self):
+        self.assertTrue(DC.Examples.Verification.run() < 0.1)
+
+
+
+if __name__ == '__main__':
+    unittest.main()