diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py
index b2a75ed4..3b5d8333 100644
--- a/SimPEG/DiffOperators.py
+++ b/SimPEG/DiffOperators.py
@@ -1,6 +1,6 @@
 import numpy as np
 from scipy import sparse as sp
-from utils import mkvc, sdiag, speye, kron3, spzeros
+from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros
 
 
 def ddx(n):
@@ -287,15 +287,19 @@ class DiffOperators(object):
     nodalVectorAve = property(**nodalVectorAve())
 
     def getEdgeMass(self, materialProp=None):
-        """mass matix for products of edge functions w'*M(materialProp)*e"""
+        """mass matrix for products of edge functions w'*M(materialProp)*e"""
         if(materialProp is None):
             materialProp = np.ones(self.nC)
         Av = self.edgeAve
         return sdiag(Av.T * (self.vol * mkvc(materialProp)))
 
     def getFaceMass(self, materialProp=None):
-        """mass matix for products of edge functions w'*M(materialProp)*e"""
+        """mass matrix for products of face functions w'*M(materialProp)*f"""
         if(materialProp is None):
             materialProp = np.ones(self.nC)
         Av = self.faceAve
-        return sdiag(Av.T*(self.vol*mkvc(materialProp)))
+        return sdiag(Av.T * (self.vol * mkvc(materialProp)))
+
+    def getFaceMassDeriv(self):
+        Av = self.faceAve
+        return Av.T * sdiag(self.vol)
diff --git a/SimPEG/GaussNewton.py b/SimPEG/GaussNewton.py
deleted file mode 100644
index e0304a8d..00000000
--- a/SimPEG/GaussNewton.py
+++ /dev/null
@@ -1,141 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from pylab import norm
-
-def GaussNewton(fctn, x0,maxIter=20, maxIterLS=10, LSreduction=1e-4, tolJ=1e-3, tolX=1e-3, 
-                                                    tolG=1e-3, eps=1e-16, xStop=[]):
-    """
-        GaussNewton Optimization
-        
-        Input:
-        ------
-            fctn - objective Function (lambda function)
-            x0   - starting guess
-            
-        Output:
-        -------
-            xOpt - numerical optimizer
-    """
-    # initial output
-    print "%s GaussNewton %s" % ('='*22,'='*22)
-    print "iter\tJc\t\tnorm(dJ)\tLS" 
-    print "%s" % '-'*57
-    
-    # evaluate stopping criteria
-    if xStop==[]:
-        xStop=x0
-    Jstop = fctn(xStop)
-    print "%3d\t%1.2e" % (-1, Jstop[0])
-
-    # initialize
-    xc = x0
-    STOP = np.zeros((5,1),dtype=bool)
-    iterLS=0; iter=0
-    
-    Jold = Jstop
-    xOld=xc
-    while 1:
-        # evaluate objective function
-        Jc,dJ,H = fctn(xc)
-        print "%3d\t%1.2e\t%1.2e\t%d" % (iter, Jc[0],norm(dJ),iterLS)
-        
-        # check stopping rules
-        STOP[0] = (iter>0) & (abs(Jc[0]-Jold[0]) <= tolJ*(1+abs(Jstop[0])))
-        STOP[1] = (iter>0) & (norm(xc-xOld)      <= tolX*(1+norm(x0)))
-        STOP[2] = norm(dJ)                       <= tolG*(1+abs(Jstop[0]))
-        STOP[3] = norm(dJ)                       <= 1e3*eps
-        STOP[4] = (iter >= maxIter)
-        if all(STOP[0:3]) | any(STOP[3:]): 
-            break
-        
-        # get search direction
-        dx = np.linalg.solve(H,-dJ)
-        
-        # Armijo linesearch
-        descent = np.dot(dJ.T,dx)
-        LS =0; t = 1; iterLS=1
-        while  (iterLS<maxIterLS):
-            xt = xc + t*dx
-            Jt = fctn(xt)
-            LS = Jt[0]<Jc[0]+t*LSreduction*descent
-            if LS:
-                break
-            iterLS = iterLS+1
-            t = .5*t
-            
-        # store old values
-        Jold = Jc; xOld = xc
-        # update 
-        xc = xt
-        iter = iter +1
-        
-    print "%s STOP! %s" % ('-'*25,'-'*25)
-    print "%d : |Jc-Jold| = %1.4e <= tolJ*(1+|Jstop|) = %1.4e"  % (STOP[0],abs(Jc[0]-Jold[0]),tolJ*(1+abs(Jstop[0])))
-    print "%d : |xc-xOld| = %1.4e <= tolX*(1+|x0|)    = %1.4e"  % (STOP[1],norm(xc-xOld),tolX*(1+norm(x0)))
-    print "%d : |dJ|      = %1.4e <= tolG*(1+|Jstop|) = %1.4e"  % (STOP[2],norm(dJ),tolG*(1+abs(Jstop[0])))
-    print "%d : |dJ|      = %1.4e <= 1e3*eps          = %1.4e"  % (STOP[3],norm(dJ),1e3*eps)
-    print "%d : iter      = %3d\t <= maxIter\t       = %3d"     % (STOP[4],iter,maxIter)
-    print "%s DONE! %s\n" % ('='*25,'='*25)
-    
-    return xc
-      
-def Rosenbrock(x):
-    """
-        Rosenbrock function for testing GaussNewton scheme
-    """
-    J   = 100*(x[1]-x[0]**2)**2+(1-x[0])**2 
-    dJ  = np.array([2*(200*x[0]**3-200*x[0]*x[1]+x[0]-1),200*(x[1]-x[0]**2)])
-    H = np.array([[-400*x[1]+1200*x[0]**2+2, -400*x[0]],[ -400*x[0], 200]],dtype=float);
-    
-    return J,dJ,H
-    
-def checkDerivative(fctn,x0):
-    """
-        Basic derivative check
-        
-        Compares error decay of 0th and 1st order Taylor approximation at point 
-        x0 for a randomized search direction.
-        
-       Input:
-       ------
-         fctn  -  function handle
-         x0    -  point at which to check derivative   
-    """
-    
-    print "%s checkDerivative %s" % ('='*20,'='*20)
-    print "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|"
-
-    Jc,dJ,H = fctn(x0)
-    
-    dx = np.random.randn(len(x0),1)
-    
-    t  = np.logspace(-1,-10,10)
-    E0 = np.zeros(t.shape)
-    E1 = np.zeros(t.shape)
-    
-    for i in range(0,10):
-        Jt = fctn(x0+t[i]*dx)
-        E0[i] = norm(Jt[0]-Jc[0])                           # 0th order Taylor 
-        E1[i] = norm(Jt[0]-Jc[0]-t[i]*np.dot(dJ.T,dx))      # 1st order Taylor 
-        
-        print "%d\t%1.2e\t%1.3e\t%1.3e" % (i,t[i],E0[i],E1[i])
-        
-    print "%s DONE! %s\n" % ('='*25,'='*25)
-    plt.figure()
-    plt.clf()
-    plt.loglog(t,E0,'b')
-    plt.loglog(t,E1,'g--')
-    plt.title('checkDerivative')
-    plt.xlabel('h')
-    plt.ylabel('error of Taylor approximation')
-    plt.legend(['0th order', '1st order'],loc='upper left')
-    plt.show()
-    return
-    
-if __name__ == '__main__':
-    x0 = np.array([[2.6],[3.7]])
-    fctn = lambda x:Rosenbrock(x)
-    checkDerivative(fctn,x0)
-    xOpt = GaussNewton(fctn,x0,maxIter=20)
-    print "xOpt=[%f,%f]" % (xOpt[0],xOpt[1])
-   
\ No newline at end of file
diff --git a/SimPEG/InnerProducts.py b/SimPEG/InnerProducts.py
index fca632f4..9fe84ac3 100644
--- a/SimPEG/InnerProducts.py
+++ b/SimPEG/InnerProducts.py
@@ -1,5 +1,5 @@
 from scipy import sparse as sp
-from utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
+from SimPEG.utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
 import numpy as np
 
 
diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py
new file mode 100644
index 00000000..2087ef5f
--- /dev/null
+++ b/SimPEG/Solver.py
@@ -0,0 +1,75 @@
+import numpy as np
+import scipy.sparse.linalg as linalg
+
+
+class Solver(object):
+    """docstring for Solver"""
+    def __init__(self, A, doDirect=True, flag=None, options={}):
+
+        assert type(doDirect) is bool, 'doDirect must be a boolean'
+        assert flag in [None, 'L', 'U', 'D'], "flag must be set to None, 'L', 'U', or 'D'"
+
+        self.A = A
+
+        self.dsolve = None
+        self.doDirect = doDirect
+        self.flag = flag
+        self.options = options
+
+    def solve(self, b):
+        if self.flag is None and self.doDirect:
+            return self.solveDirect(b, **self.options)
+        elif self.flag is None and not self.doDirect:
+            return self.solveIter(b, **self.options)
+        elif self.flag == 'U':
+            return self.solveBackward(b)
+        elif self.flag == 'L':
+            return self.solveForward(b)
+        elif self.flag == 'D':
+            return self.solveDiagonal(b)
+        else:
+            raise Exception('Unknown flag.')
+        pass
+
+    def clean(self):
+        """Cleans up the memory"""
+        del self.dsolve
+        self.dsolve = None
+
+    def solveDirect(self, b, backend='scipy'):
+        assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch'
+
+        if self.dsolve is None:
+            self.A = self.A.tocsc()  # for efficiency
+            self.dsolve = linalg.factorized(self.A)
+
+        if len(b.shape) == 1 or b.shape[1] == 1:
+            # Just one RHS
+            return self.dsolve(b)
+
+        # Multiple RHSs
+        X = np.empty_like(b)
+        for i in range(b.shape[1]):
+            X[:,i] = self.dsolve(b[:,i])
+
+        return X
+
+    def solveIter(self, b, M=None, iterSolver='CG'):
+        pass
+
+    def solveBackward(self, b):
+        pass
+
+    def solveForward(self, b):
+        pass
+
+    def solveDiagonal(self, b):
+        diagA = self.A.diagonal()
+        if len(b.shape) == 1 or b.shape[1] == 1:
+            # Just one RHS
+            return b/diagA
+        # Multiple RHSs
+        X = np.empty_like(b)
+        for i in range(b.shape[1]):
+            X[:,i] = b[:,i]/diagA
+        return X
diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py
index ca9b9f7b..ea9faf33 100644
--- a/SimPEG/__init__.py
+++ b/SimPEG/__init__.py
@@ -1,3 +1,5 @@
 from TensorMesh import TensorMesh
 from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh
 import utils
+import inverse
+from Solver import Solver
diff --git a/SimPEG/forward/DCProblem/DCProblem.py b/SimPEG/forward/DCProblem/DCProblem.py
new file mode 100644
index 00000000..b8115598
--- /dev/null
+++ b/SimPEG/forward/DCProblem/DCProblem.py
@@ -0,0 +1,168 @@
+from SimPEG import TensorMesh
+from SimPEG.forward import Problem, SyntheticProblem
+from SimPEG.tests import checkDerivative
+from SimPEG.utils import ModelBuilder, sdiag
+import numpy as np
+import scipy.sparse.linalg as linalg
+import DCutils
+
+class DCProblem(Problem):
+    """
+        **DCProblem**
+
+        Geophysical DC resistivity problem.
+
+    """
+    def __init__(self, mesh):
+        super(DCProblem, self).__init__(mesh)
+        self.mesh.setCellGradBC('neumann')
+
+    def createMatrix(self, m):
+        """
+            Makes the matrix A(m) for the DC resistivity problem.
+
+            :param numpy.array m: model
+            :rtype: scipy.csc_matrix
+            :return: A(m)
+
+            .. math::
+                c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
+
+            Where M() is the mass matrix and mT is the model transform.
+        """
+        D = self.mesh.faceDiv
+        G = self.mesh.cellGrad
+        sigma = self.modelTransform(m)
+        Msig = self.mesh.getFaceMass(sigma)
+        A = D*Msig*G
+        return A.tocsc()
+
+    def field(self, m):
+        A = self.createMatrix(m)
+        solve = linalg.factorized(A)
+
+        nRHSs = self.RHS.shape[1]  # Number of RHSs
+        phi = np.zeros((self.mesh.nC, nRHSs)) + np.nan
+        for ii in range(nRHSs):
+            phi[:,ii] = solve(self.RHS[:,ii])
+
+        return phi
+
+    def J(self, m, v, u=None, solve=None):
+        """
+            :param numpy.array m: model
+            :param numpy.array v: vector to multiply
+            :param numpy.array u: fields
+            :rtype: numpy.array
+            :return: Jv
+
+            .. math::
+                c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
+
+                \\nabla_u (A(m)u - q) = A(m)
+
+                \\nabla_m (A(m)u - q) = G\\text{sdiag}(Du)\\nabla_m(M(mT(m)))
+
+            Where M() is the mass matrix and mT is the model transform.
+
+            .. math::
+                J = - P \left( \\nabla_u c(m, u) \\right)^{-1} \\nabla_m c(m, u)
+
+                J(v) = - P ( A(m)^{-1} ( G\\text{sdiag}(Du)\\nabla_m(M(mT(m))) v ) )
+        """
+        P = self.P
+        D = self.mesh.faceDiv
+        G = self.mesh.cellGrad
+        A = self.createMatrix(m)
+        Av_dm = self.mesh.getFaceMassDeriv()
+        mT_dm = self.modelTransformDeriv(m)
+
+        dCdu = A
+        dCdm = D * ( sdiag( G * u ) * ( Av_dm * ( mT_dm * v ) ) )
+
+        if solve is None:
+            solve = linalg.factorized(dCdu)
+
+        Jv = - P * solve(dCdm)
+        return Jv
+
+    def Jt(self, m, v, u=None, solve=None):
+        P = self.P
+        D = self.mesh.faceDiv
+        G = self.mesh.cellGrad
+        A = self.createMatrix(m)
+        Av_dm = self.mesh.getFaceMassDeriv()
+        mT_dm = self.modelTransformDeriv(m)
+
+        dCdu = A.T
+
+        if solve is None:
+            solve = linalg.factorized(dCdu.tocsc())
+        w = solve(P.T*v)
+
+        Jtv = - mT_dm.T * ( Av_dm.T * ( sdiag( G * u ) * ( D.T * w ) ) )
+        return Jtv
+
+
+if __name__ == '__main__':
+    # Create the mesh
+    h1 = np.ones(100)
+    h2 = np.ones(100)
+    mesh = TensorMesh([h1,h2])
+
+    # Create some parameters for the model
+    sig1 = 1
+    sig2 = 0.01
+
+    # Create a synthetic model from a block in a half-space
+    p0 = [20, 20]
+    p1 = [50, 50]
+    condVals = [sig1, sig2]
+    mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals)
+    mesh.plotImage(mSynth, showIt=False)
+
+
+    # Set up the projection
+    nelec = 50
+    spacelec = 2
+    surfloc = 0.5
+    elecini = 0.5
+    elecend = 0.5+spacelec*(nelec-1)
+    elecLocR = np.linspace(elecini, elecend, nelec)
+    rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
+    q, Q, rxmidloc = DCutils.genTxRxmat(nelec, spacelec, surfloc, elecini, mesh)
+    P = Q.T
+
+    # Create some data
+    class syntheticDCProblem(DCProblem, SyntheticProblem):
+        pass
+
+    synthetic = syntheticDCProblem(mesh);
+    synthetic.P = P
+    synthetic.RHS = q
+    dobs, Wd = synthetic.createData(mSynth, std=0.05)
+
+    u = synthetic.field(mSynth)
+    mesh.plotImage(u[:,10], showIt=True)
+
+    # Now set up the problem to do some minimization
+    problem = DCProblem(mesh)
+    problem.P = P
+    problem.RHS = q
+    problem.W = Wd
+    problem.dobs = dobs
+    m0 = mesh.gridCC[:,0]*0+sig1
+
+    print problem.misfit(m0)
+    print problem.misfit(mSynth)
+
+    # Check Derivative
+    derChk = lambda m: [problem.misfit(m), problem.misfitDeriv(m)]
+    checkDerivative(derChk, mSynth)
+
+    # Adjoint Test
+    u = np.random.rand(mesh.nC)
+    v = np.random.rand(mesh.nC)
+    w = np.random.rand(dobs.shape[0])
+    print w.dot(problem.J(mSynth, v, u=u))
+    print v.dot(problem.Jt(mSynth, w, u=u))
diff --git a/SimPEG/forward/DCProblem/DCutils.py b/SimPEG/forward/DCProblem/DCutils.py
new file mode 100644
index 00000000..f3445096
--- /dev/null
+++ b/SimPEG/forward/DCProblem/DCutils.py
@@ -0,0 +1,29 @@
+import numpy as np
+import scipy.sparse as sp
+
+def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh):
+    """ Generate projection matrix (Q) and """
+    elecend = 0.5+spacelec*(nelec-1)
+    elecLocR = np.linspace(elecini, elecend, nelec)
+    elecLocT = elecLocR+1
+    nrx = nelec-1
+    ntx = nelec-1
+    q = np.zeros((mesh.nC, ntx))
+    Q = np.zeros((mesh.nC, nrx))
+
+    for i in range(nrx):
+
+        rxind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i]))
+        rxind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i+1]))
+
+        txind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i]))
+        txind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i+1]))
+
+        q[txind1,i] = 1
+        q[txind2,i] = -1
+        Q[rxind1,i] = 1
+        Q[rxind2,i] = -1
+
+    Q = sp.csr_matrix(Q)
+    rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
+    return q, Q, rxmidLoc
diff --git a/SimPEG/forward/DCProblem/__init__.py b/SimPEG/forward/DCProblem/__init__.py
new file mode 100644
index 00000000..a868cf80
--- /dev/null
+++ b/SimPEG/forward/DCProblem/__init__.py
@@ -0,0 +1,2 @@
+from DCProblem import *
+from DCutils import *
diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py
new file mode 100644
index 00000000..5b716f1f
--- /dev/null
+++ b/SimPEG/forward/Problem.py
@@ -0,0 +1,339 @@
+import numpy as np
+from SimPEG.utils import mkvc, sdiag
+norm = np.linalg.norm
+
+
+class Problem(object):
+    """
+        Problem is the base class for all geophysical forward problems in SimPEG.
+
+
+        The problem is a partial differential equation of the form:
+
+        .. math::
+            c(m, u) = 0
+
+        Here, m is the model and u is the field (or fields).
+        Given the model, m, we can calculate the fields u(m),
+        however, the data we collect is a subset of the fields,
+        and can be defined by a linear projection, P.
+
+        .. math::
+            d_\\text{pred} = Pu(m)
+
+        We are interested in how changing the model transforms the data,
+        as such we can take write the Taylor expansion:
+
+        .. math::
+            Pu(m + hv) = Pu(m) + hP\\frac{\partial u(m)}{\partial m} v + \mathcal{O}(h^2 \left\| v \\right\| )
+
+        We can linearize and define the sensitivity matrix as:
+
+        .. math::
+            J = P\\frac{\partial u}{\partial m}
+
+        The sensitivity matrix, and it's transpose will be used in the inverse problem
+        to (locally) find how model parameters change the data, and optimize!
+    """
+
+    def __init__(self, mesh):
+        self.mesh = mesh
+
+    @property
+    def RHS(self):
+        """
+            Source matrix.
+        """
+        return self._RHS
+    @RHS.setter
+    def RHS(self, value):
+        self._RHS = value
+
+    @property
+    def W(self):
+        """
+            Standard deviation weighting matrix.
+        """
+        return self._W
+    @W.setter
+    def W(self, value):
+        self._W = value
+
+    @property
+    def P(self):
+        """
+            Projection matrix.
+
+            .. math::
+                d_\\text{pred} = Pu(m)
+        """
+        return self._P
+    @P.setter
+    def P(self, value):
+        self._P = value
+
+
+    @property
+    def dobs(self):
+        """
+            Observed data.
+        """
+        return self._dobs
+    @dobs.setter
+    def dobs(self, value):
+        self._dobs = value
+
+    def evalFunction(self, m, doDerivative=True):
+        """
+            :param numpy.array m: model
+            :param bool doDerivative: do you want to compute the derivative?
+            :rtype: numpy.array
+            :return: Jv
+        """
+        f = self.misfit(m)
+
+        return f, g, H
+
+    def J(self, m, v, u=None):
+        """
+            :param numpy.array m: model
+            :param numpy.array v: vector to multiply
+            :param numpy.array u: fields
+            :rtype: numpy.array
+            :return: Jv
+
+            Working with the general PDE, c(m, u) = 0, where m is the model and u is the field,
+            the sensitivity is defined as:
+
+            .. math::
+                J = P\\frac{\partial u}{\partial m}
+
+            We can take the derivative of the PDE:
+
+            .. math::
+                \\nabla_m c(m, u) \delta m + \\nabla_u c(m, u) \delta u = 0
+
+            If the forward problem is invertible, then we can rearrange for du/dm:
+
+            .. math::
+                J = - P \left( \\nabla_u c(m, u) \\right)^{-1} \\nabla_m c(m, u)
+
+            This can often be computed given a vector (i.e. J(v)) rather than stored, as J is a large dense matrix.
+
+        """
+        pass
+
+    def Jt(self, m, v, u=None):
+        """
+            :param numpy.array m: model
+            :param numpy.array v: vector to multiply
+            :param numpy.array u: fields
+            :rtype: numpy.array
+            :return: JTv
+
+            Transpose of J
+        """
+        pass
+
+    def field(self, m):
+        """
+            The field given the model.
+
+            .. math::
+                u(m)
+
+        """
+        pass
+
+    def dpred(self, m, u=None):
+        """
+            Predicted data.
+
+            .. math::
+                d_\\text{pred} = Pu(m)
+        """
+        if u is None:
+            u = self.field(m)
+        return self.P*u
+
+    def modelTransform(self, m):
+        """
+            :param numpy.array m: model
+            :rtype: numpy.array
+            :return: transformed model
+
+            The modelTransform changes the model into the physical property.
+
+            A common example of this is to invert for electrical conductivity
+            in log space. In this case, your model will be log(sigma) and to
+            get back to sigma, you can take the exponential:
+
+            .. math::
+
+                m = \log{\sigma}
+
+                \exp{m} = \exp{\log{\sigma}} = \sigma
+        """
+        return np.exp(mkvc(m))
+
+    def modelTransformDeriv(self, m):
+        """
+            :param numpy.array m: model
+            :rtype: scipy.csr_matrix
+            :return: derivative of transformed model
+
+            The modelTransform changes the model into the physical property.
+            The modelTransformDeriv provides the derivative of the modelTransform.
+
+            If the model transform is:
+
+            .. math::
+
+                m = \log{\sigma}
+
+                \exp{m} = \exp{\log{\sigma}} = \sigma
+
+            Then the derivative is:
+
+            .. math::
+
+                \\frac{\partial \exp{m}}{\partial m} = \\text{sdiag}(\exp{m})
+        """
+        return sdiag(np.exp(mkvc(m)))
+
+    def misfit(self, m, u=None):
+        """
+            :param numpy.array m: geophysical model
+            :param numpy.array u: fields
+            :rtype: float
+            :return: data misfit
+
+            The data misfit using an l_2 norm is:
+
+            .. math::
+
+                \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2
+
+            Where P is a projection matrix that brings the field on the full domain to the data measurement locations;
+            u is the field of interest; d_obs is the observed data; and W is the weighting matrix.
+        """
+
+        R = self.W*(self.dpred(m, u=u) - self.dobs)
+        R = mkvc(R)
+        return 0.5*R.dot(R)
+
+    def misfitDeriv(self, m, u=None):
+        """
+            :param numpy.array m: geophysical model
+            :param numpy.array u: fields
+            :rtype: numpy.array
+            :return: data misfit derivative
+
+            The data misfit using an l_2 norm is:
+
+            .. math::
+
+                \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2
+
+            If the field, u, is provided, the calculation of the data is fast:
+
+            .. math::
+
+                \mathbf{d}_\\text{pred} = \mathbf{Pu(m)}
+
+                \mathbf{R} = \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs})
+
+            Where P is a projection matrix that brings the field on the full domain to the data measurement locations;
+            u is the field of interest; d_obs is the observed data; and W is the weighting matrix.
+
+            The derivative of this, with respect to the model, is:
+
+            .. math::
+
+                \\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R}
+
+        """
+        if u is None:
+            u = self.field(m)
+
+        R = self.W*(self.dpred(m, u=u) - self.dobs)
+
+        dmisfit = 0
+        for i in range(self.RHS.shape[1]): # Loop over each right hand side
+            dmisfit += self.Jt(m, self.W[:,i]*R[:,i], u=u[:,i])
+
+        return dmisfit
+
+    def misfitDerivDeriv(self, m, u=None):
+        """
+            :param numpy.array m: geophysical model
+            :param numpy.array u: fields
+            :rtype: numpy.array
+            :return: data misfit derivative
+
+            The data misfit using an l_2 norm is:
+
+            .. math::
+
+                \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2
+
+            If the field, u, is provided, the calculation of the data is fast:
+
+            .. math::
+
+                \mathbf{d}_\\text{pred} = \mathbf{Pu(m)}
+
+                \mathbf{R} = \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs})
+
+            Where P is a projection matrix that brings the field on the full domain to the data measurement locations;
+            u is the field of interest; d_obs is the observed data; and W is the weighting matrix.
+
+            The derivative of this, with respect to the model, is:
+
+            .. math::
+
+                \\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R}
+
+                \\frac{\partial^2 \mu_\\text{data}}{\partial^2 \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ W J}
+
+        """
+        if u is None:
+            u = self.field(m)
+
+        R = self.W*(self.dpred(m, u=u) - self.dobs)
+
+        dmisfit = 0
+        for i in range(self.RHS.shape[1]): # Loop over each right hand side
+            dmisfit += self.Jt(m, self.W[:,i]*R[:,i], u=u[:,i])
+
+        return dmisfit
+
+
+class SyntheticProblem(object):
+    """
+        Has helpful functions when dealing with synthetic problems
+
+        To use this class, inherit to your problem::
+
+            class mySyntheticExample(Problem, SyntheticProblem):
+                pass
+    """
+    def createData(self, m, std=0.05):
+        """
+            :param numpy.array m: geophysical model
+            :param numpy.array std: standard deviation
+            :rtype: numpy.array, numpy.array
+            :return: dobs, Wd
+
+            Create synthetic data given a model, and a standard deviation.
+
+            Returns the observed data with random Gaussian noise
+            and Wd which is the same size as data, and can be used to weight the inversion.
+        """
+        dobs = self.dpred(m)
+        dobs = dobs
+        noise = std*abs(dobs)*np.random.randn(*dobs.shape)
+        dobs = dobs+noise
+        eps = np.linalg.norm(mkvc(dobs),2)*1e-5
+        Wd = 1/(abs(dobs)*std+eps)
+        return dobs, Wd
diff --git a/SimPEG/forward/__init__.py b/SimPEG/forward/__init__.py
new file mode 100644
index 00000000..fe849d41
--- /dev/null
+++ b/SimPEG/forward/__init__.py
@@ -0,0 +1,2 @@
+from Problem import *
+import DCProblem
diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py
new file mode 100644
index 00000000..eee18b16
--- /dev/null
+++ b/SimPEG/inverse/Optimize.py
@@ -0,0 +1,150 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from SimPEG.utils import mkvc, sdiag
+norm = np.linalg.norm
+
+
+class Minimize(object):
+    """docstring for Minimize"""
+
+    name = "GeneralOptimizationAlgorithm"
+
+    maxIter = 20
+    maxIterLS = 10
+    LSreduction = 1e-4
+    LSshorten = 0.5
+    tolF = 1e-4
+    tolX = 1e-4
+    tolG = 1e-4
+    eps = 1e-16
+
+    def __init__(self, problem, **kwargs):
+        self.problem = problem
+        self.setKwargs(**kwargs)
+
+    def setKwargs(self, **kwargs):
+        # Set the variables, throw an error if they don't exist.
+        for attr in kwargs:
+            if hasattr(self, attr):
+                setattr(self, attr, kwargs[attr])
+            else:
+                raise Exception('%s attr is not recognized' % attr)
+
+    def minimize(self, x0):
+
+        self.startup(x0)
+        self.printInit()
+
+        while True:
+            self.f, self.g, self.H = self.evalFunction(self.xc)
+            self.printIter()
+            if self.stoppingCriteria(): break
+            p = self.findSearchDirection()
+            xt, passLS = self.linesearch(p)
+            if not passLS:
+                xt = self.linesearchBreak(p)
+            self.doEndIteration(xt)
+
+        self.printDone()
+
+        return self.xc
+
+    def startup(self, x0):
+        self._iter = 0
+        self._iterLS = 0
+        self._STOP = np.zeros((5,1),dtype=bool)
+
+        self.x0 = x0
+        self.xc = x0
+        self.xOld = x0
+
+    def printInit(self):
+        print "%s %s %s" % ('='*22, self.name, '='*22)
+        print "iter\tJc\t\tnorm(dJ)\tLS"
+        print "%s" % '-'*57
+
+    def printIter(self):
+        print "%3d\t%1.2e\t%1.2e\t%d" % (self._iter, self.f, norm(self.g), self._iterLS)
+
+    def printDone(self):
+        print "%s STOP! %s" % ('-'*25,'-'*25)
+        print "%d : |fc-fOld| = %1.4e <= tolF*(1+|fStop|) = %1.4e"  % (self._STOP[0], abs(self.f-self.fOld), self.tolF*(1+abs(self.fStop)))
+        print "%d : |xc-xOld| = %1.4e <= tolX*(1+|x0|)    = %1.4e"  % (self._STOP[1], norm(self.xc-self.xOld), self.tolX*(1+norm(self.x0)))
+        print "%d : |g|       = %1.4e <= tolG*(1+|fStop|) = %1.4e"  % (self._STOP[2], norm(self.g), self.tolG*(1+abs(self.fStop)))
+        print "%d : |g|       = %1.4e <= 1e3*eps          = %1.4e"  % (self._STOP[3], norm(self.g), 1e3*self.eps)
+        print "%d : iter      = %3d\t <= maxIter\t       = %3d"     % (self._STOP[4], self._iter, self.maxIter)
+        print "%s DONE! %s\n" % ('='*25,'='*25)
+
+    def evalFunction(self, x, doDerivative=True):
+        f, g, H = self.problem(x)
+        return f, g, H
+
+    def findSearchDirection(self):
+        return -self.g
+
+    def stoppingCriteria(self):
+        if self._iter == 0:
+            self.fStop = self.f  # Save this for stopping criteria
+
+        # check stopping rules
+        self._STOP[0] = self._iter > 0 and (abs(self.f-self.fOld) <= self.tolF*(1+abs(self.fStop)))
+        self._STOP[1] = self._iter > 0 and (norm(self.xc-self.xOld) <= self.tolX*(1+norm(self.x0)))
+        self._STOP[2] = norm(self.g) <= self.tolG*(1+abs(self.fStop))
+        self._STOP[3] = norm(self.g) <= 1e3*self.eps
+        self._STOP[4] = self._iter >= self.maxIter
+        return all(self._STOP[0:3]) | any(self._STOP[3:])
+
+    def linesearch(self, p):
+        # Armijo linesearch
+        descent = np.inner(self.g, p)
+        t = 1
+        iterLS = 0
+        while iterLS < self.maxIterLS:
+            xt = self.xc + t*p
+            ft, temp, temp = self.evalFunction(xt, doDerivative=False)
+            if ft < self.f + t*self.LSreduction*descent:
+                break
+            iterLS += 1
+            t = self.LSshorten*t
+
+        self._iterLS = iterLS
+        return xt, iterLS < self.maxIterLS
+
+    def linesearchBreak(self, p):
+        raise Exception('The linesearch got broken. Boo.')
+
+    def doEndIteration(self, xt):
+        # store old values
+        self.fOld = self.f
+        self.xOld, self.xc = self.xc, xt
+        self._iter += 1
+
+
+class GaussNewton(Minimize):
+    name = 'GaussNewton'
+    def findSearchDirection(self):
+        return np.linalg.solve(self.H,-self.g)
+
+
+class SteepestDescent(Minimize):
+    name = 'SteepestDescent'
+    def findSearchDirection(self):
+        return -self.g
+
+if __name__ == '__main__':
+    from SimPEG.tests import Rosenbrock, checkDerivative
+    x0 = np.array([2.6, 3.7])
+    checkDerivative(Rosenbrock, x0, plotIt=False)
+    xOpt = GaussNewton(Rosenbrock, maxIter=20).minimize(x0)
+    print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1])
+    xOpt = SteepestDescent(Rosenbrock, maxIter=20, maxIterLS=15).minimize(x0)
+    print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1])
+
+    def simplePass(x):
+        return np.sin(x), sdiag(np.cos(x))
+
+    def simpleFail(x):
+        return np.sin(x), -sdiag(np.cos(x))
+
+    checkDerivative(simplePass, np.random.randn(5), plotIt=False)
+    checkDerivative(simpleFail, np.random.randn(5), plotIt=False)
diff --git a/SimPEG/inverse/__init__.py b/SimPEG/inverse/__init__.py
new file mode 100644
index 00000000..b2a5e506
--- /dev/null
+++ b/SimPEG/inverse/__init__.py
@@ -0,0 +1 @@
+from Optimize import *
diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/TestUtils.py
similarity index 56%
rename from SimPEG/tests/OrderTest.py
rename to SimPEG/tests/TestUtils.py
index fb305d32..79340da7 100644
--- a/SimPEG/tests/OrderTest.py
+++ b/SimPEG/tests/TestUtils.py
@@ -1,5 +1,7 @@
-import sys
-sys.path.append('../../')
+import numpy as np
+import matplotlib.pyplot as plt
+from pylab import norm
+from SimPEG.utils import mkvc
 from SimPEG import TensorMesh, utils, LogicallyOrthogonalMesh
 import numpy as np
 import unittest
@@ -16,46 +18,47 @@ class OrderTest(unittest.TestCase):
         Given are an operator A and its discretization A[h]. For a given test function f
         and h --> 0  we compare:
 
-        error(h) = \| A[h](f) - A(f) \|_{\infty}
+        .. math::
+            error(h) = \| A[h](f) - A(f) \|_{\infty}
 
         Note that you can provide any norm.
 
         Test is passed when estimated rate order of convergence is  at least within the specified tolerance of the
         estimated rate supplied by the user.
 
-    Minimal example for a curl operator:
+    Minimal example for a curl operator::
 
-    class TestCURL(OrderTest):
-        name = "Curl"
+        class TestCURL(OrderTest):
+            name = "Curl"
 
-        def getError(self):
-            # For given Mesh, generate A[h], f and A(f) and return norm of error.
+            def getError(self):
+                # For given Mesh, generate A[h], f and A(f) and return norm of error.
 
 
-            fun  = lambda x: np.cos(x)  # i (cos(y)) + j (cos(z)) + k (cos(x))
-            sol = lambda x: np.sin(x)  # i (sin(z)) + j (sin(x)) + k (sin(y))
+                fun  = lambda x: np.cos(x)  # i (cos(y)) + j (cos(z)) + k (cos(x))
+                sol = lambda x: np.sin(x)  # i (sin(z)) + j (sin(x)) + k (sin(y))
 
 
-            Ex = fun(self.M.gridEx[:, 1])
-            Ey = fun(self.M.gridEy[:, 2])
-            Ez = fun(self.M.gridEz[:, 0])
-            f = np.concatenate((Ex, Ey, Ez))
+                Ex = fun(self.M.gridEx[:, 1])
+                Ey = fun(self.M.gridEy[:, 2])
+                Ez = fun(self.M.gridEz[:, 0])
+                f = np.concatenate((Ex, Ey, Ez))
 
-            Fx = sol(self.M.gridFx[:, 2])
-            Fy = sol(self.M.gridFy[:, 0])
-            Fz = sol(self.M.gridFz[:, 1])
-            Af = np.concatenate((Fx, Fy, Fz))
+                Fx = sol(self.M.gridFx[:, 2])
+                Fy = sol(self.M.gridFy[:, 0])
+                Fz = sol(self.M.gridFz[:, 1])
+                Af = np.concatenate((Fx, Fy, Fz))
 
-            # Generate DIV matrix
-            Ah = self.M.edgeCurl
+                # Generate DIV matrix
+                Ah = self.M.edgeCurl
 
-            curlE = Ah*E
-            err = np.linalg.norm((Ah*f -Af), np.inf)
-            return err
+                curlE = Ah*E
+                err = np.linalg.norm((Ah*f -Af), np.inf)
+                return err
 
-        def test_order(self):
-            # runs the test
-            self.orderTest()
+            def test_order(self):
+                # runs the test
+                self.orderTest()
 
     See also: test_operatorOrder.py
 
@@ -159,5 +162,78 @@ class OrderTest(unittest.TestCase):
             print ''
             self.assertTrue(passTest)
 
-if __name__ == '__main__':
-    unittest.main()
+def Rosenbrock(x):
+    """Rosenbrock function for testing GaussNewton scheme"""
+
+    f = 100*(x[1]-x[0]**2)**2+(1-x[0])**2
+    g = np.array([2*(200*x[0]**3-200*x[0]*x[1]+x[0]-1), 200*(x[1]-x[0]**2)])
+    H = np.array([[-400*x[1]+1200*x[0]**2+2, -400*x[0]], [-400*x[0], 200]])
+    return f, g, H
+
+def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None):
+    """
+        Basic derivative check
+
+        Compares error decay of 0th and 1st order Taylor approximation at point
+        x0 for a randomized search direction.
+
+        :param lambda fctn: function handle
+        :param numpy.array x0: point at which to check derivative
+        :param int num: number of times to reduce step length, h
+        :param bool plotIt: if you would like to plot
+        :param numpy.array dx: step direction
+        :rtype: bool
+        :return: did you pass the test?!
+
+    """
+
+    print "%s checkDerivative %s" % ('='*20, '='*20)
+    print "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|\tOrder\n%s" % ('-'*57)
+
+    Jc = fctn(x0)
+
+    x0 = mkvc(x0)
+
+    if dx is None:
+        dx = np.random.randn(len(x0))
+
+    t  = np.logspace(-1, -num, num)
+    E0 = np.ones(t.shape)
+    E1 = np.ones(t.shape)
+
+    l2norm = lambda x: np.sqrt(np.inner(x, x))  # because np.norm breaks if they are scalars?
+    for i in range(num):
+        Jt = fctn(x0+t[i]*dx)
+        E0[i] = l2norm(Jt[0]-Jc[0])               # 0th order Taylor
+        E1[i] = l2norm(Jt[0]-Jc[0]-t[i]*Jc[1].dot(dx))  # 1st order Taylor
+        order0 = np.log10(E0[:-1]/E0[1:])
+        order1 = np.log10(E1[:-1]/E1[1:])
+        print "%d\t%1.2e\t%1.3e\t\t%1.3e\t\t%1.3f" % (i, t[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1])
+
+    tolerance = 0.9
+    expectedOrder = 2
+    eps = 1e-10
+    order0 = order0[E0[1:] > eps]
+    order1 = order1[E1[1:] > eps]
+    belowTol = order1.size == 0 and order0.size > 0
+    correctOrder = order1.size > 0 and np.mean(order1) > tolerance * expectedOrder
+
+    passTest = belowTol or correctOrder
+
+    if passTest:
+        print "%s PASS! %s\n" % ('='*25, '='*25)
+    else:
+        print "%s\n%s FAIL! %s\n%s" % ('*'*57, '<'*25, '>'*25, '*'*57)
+
+    if plotIt:
+        plt.figure()
+        plt.clf()
+        plt.loglog(t, E0, 'b')
+        plt.loglog(t, E1, 'g--')
+        plt.title('checkDerivative')
+        plt.xlabel('h')
+        plt.ylabel('error of Taylor approximation')
+        plt.legend(['0th order', '1st order'], loc='upper left')
+        plt.show()
+
+    return passTest
diff --git a/SimPEG/tests/__init__.py b/SimPEG/tests/__init__.py
new file mode 100644
index 00000000..f896cfd0
--- /dev/null
+++ b/SimPEG/tests/__init__.py
@@ -0,0 +1,2 @@
+import TestUtils
+from TestUtils import checkDerivative, Rosenbrock, OrderTest
diff --git a/SimPEG/tests/test_forward_DCproblem.py b/SimPEG/tests/test_forward_DCproblem.py
new file mode 100644
index 00000000..b5e6d844
--- /dev/null
+++ b/SimPEG/tests/test_forward_DCproblem.py
@@ -0,0 +1,81 @@
+import numpy as np
+import unittest
+from SimPEG import TensorMesh
+from SimPEG.utils import ModelBuilder, sdiag
+from SimPEG.forward import Problem, SyntheticProblem
+from SimPEG.forward.DCProblem import DCProblem, DCutils
+from TestUtils import checkDerivative
+from scipy.sparse.linalg import dsolve
+
+
+class DCProblemTests(unittest.TestCase):
+
+    def setUp(self):
+        # Create the mesh
+        h1 = np.ones(20)
+        h2 = np.ones(20)
+        mesh = TensorMesh([h1,h2])
+
+        # Create some parameters for the model
+        sig1 = 1
+        sig2 = 0.01
+
+        # Create a synthetic model from a block in a half-space
+        p0 = [2, 2]
+        p1 = [5, 5]
+        condVals = [sig1, sig2]
+        mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals)
+
+        # Set up the projection
+        nelec = 10
+        spacelec = 2
+        surfloc = 0.5
+        elecini = 0.5
+        elecend = 0.5+spacelec*(nelec-1)
+        elecLocR = np.linspace(elecini, elecend, nelec)
+        rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
+        q, Q, rxmidloc = DCutils.genTxRxmat(nelec, spacelec, surfloc, elecini, mesh)
+        P = Q.T
+
+        # Create some data
+        class syntheticDCProblem(DCProblem, SyntheticProblem):
+            pass
+
+        synthetic = syntheticDCProblem(mesh);
+        synthetic.P = P
+        synthetic.RHS = q
+        dobs, Wd = synthetic.createData(mSynth, std=0.05)
+
+        # Now set up the problem to do some minimization
+        problem = DCProblem(mesh)
+        problem.P = P
+        problem.RHS = q
+        problem.W = Wd
+        problem.dobs = dobs
+
+        self.p = problem
+        self.mesh = mesh
+        self.m0 = mSynth
+        self.dobs = dobs
+
+
+    def test_misfit(self):
+        print 'SimPEG.forward.DCProblem: Testing Misfit'
+        derChk = lambda m: [self.p.misfit(m), self.p.misfitDeriv(m)]
+        passed = checkDerivative(derChk, self.m0, plotIt=False)
+        self.assertTrue(passed)
+
+    def test_adjoint(self):
+        # Adjoint Test
+        u = np.random.rand(self.mesh.nC)
+        v = np.random.rand(self.mesh.nC)
+        w = np.random.rand(self.dobs.shape[0])
+        wtJv = w.dot(self.p.J(self.m0, v, u=u))
+        vtJtw = v.dot(self.p.Jt(self.m0, w, u=u))
+        passed = (wtJv - vtJtw) < 1e-10
+        self.assertTrue(passed)
+
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/SimPEG/tests/test_forward_problem.py b/SimPEG/tests/test_forward_problem.py
new file mode 100644
index 00000000..b7180dca
--- /dev/null
+++ b/SimPEG/tests/test_forward_problem.py
@@ -0,0 +1,28 @@
+import numpy as np
+import unittest
+from SimPEG import TensorMesh
+from SimPEG.forward import Problem
+from TestUtils import checkDerivative
+from scipy.sparse.linalg import dsolve
+
+
+class ProblemTests(unittest.TestCase):
+
+    def setUp(self):
+
+        a = np.array([1, 1, 1])
+        b = np.array([1, 2])
+        c = np.array([1, 4])
+        self.mesh2 = TensorMesh([a, b], np.array([3, 5]))
+        self.p2 = Problem(self.mesh2)
+
+
+    def test_modelTransform(self):
+        print 'SimPEG.forward.Problem: Testing Model Transform'
+        m = np.random.rand(self.mesh2.nC)
+        passed = checkDerivative(lambda m : [self.p2.modelTransform(m), self.p2.modelTransformDeriv(m)], m, plotIt=False)
+        self.assertTrue(passed)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py
index 13037156..79626ca1 100644
--- a/SimPEG/tests/test_massMatrices.py
+++ b/SimPEG/tests/test_massMatrices.py
@@ -1,6 +1,6 @@
 import numpy as np
 import unittest
-from OrderTest import OrderTest
+from TestUtils import OrderTest
 
 
 # MATLAB code:
diff --git a/SimPEG/tests/test_operators.py b/SimPEG/tests/test_operators.py
index ab377089..6906c580 100644
--- a/SimPEG/tests/test_operators.py
+++ b/SimPEG/tests/test_operators.py
@@ -1,6 +1,6 @@
 import numpy as np
 import unittest
-from OrderTest import OrderTest
+from TestUtils import OrderTest
 
 MESHTYPES = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM']
 call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
diff --git a/SimPEG/tests/test_tensorMesh.py b/SimPEG/tests/test_tensorMesh.py
index 32d47353..300e6634 100644
--- a/SimPEG/tests/test_tensorMesh.py
+++ b/SimPEG/tests/test_tensorMesh.py
@@ -1,9 +1,7 @@
 import numpy as np
 import unittest
-import sys
-sys.path.append('../')
-from TensorMesh import TensorMesh
-from OrderTest import OrderTest
+from SimPEG import TensorMesh
+from TestUtils import OrderTest
 from scipy.sparse.linalg import dsolve
 
 
diff --git a/SimPEG/tests/test_utils.py b/SimPEG/tests/test_utils.py
index 0dcdd362..f0b867ec 100644
--- a/SimPEG/tests/test_utils.py
+++ b/SimPEG/tests/test_utils.py
@@ -1,8 +1,6 @@
 import numpy as np
 import unittest
-import sys
-sys.path.append('../')
-from utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
+from SimPEG.utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
 
 
 class TestSequenceFunctions(unittest.TestCase):
diff --git a/SimPEG/utils/ModelBuilder.py b/SimPEG/utils/ModelBuilder.py
index 3b1f977f..170da1bf 100644
--- a/SimPEG/utils/ModelBuilder.py
+++ b/SimPEG/utils/ModelBuilder.py
@@ -3,22 +3,20 @@ import numpy as np
 
 def getIndecesBlock(p0,p1,ccMesh):
     """
-    Creates a vector containing the block indexes in the cell centerd mesh.
-    Returns a tuple
+        Creates a vector containing the block indexes in the cell centerd mesh.
+        Returns a tuple
 
-    The block is defined by the points
-    p0 : describe the position of the left  upper  front corner, and
-    p1 : describe the position of the right bottom back  corner.
+        The block is defined by the points
 
-    ccMesh represents the cell-centered mesh
+        p0, describe the position of the left  upper  front corner, and
 
-    The points p0 and p1 must live in the the same dimensional space as the mesh.
+        p1, describe the position of the right bottom back  corner.
+
+        ccMesh represents the cell-centered mesh
+
+        The points p0 and p1 must live in the the same dimensional space as the mesh.
     """
 
-    # Validation of the input
-    assert type(p0) == np.ndarray, "Vector must be a numpy array"
-    assert type(p1) == np.ndarray, "Vector must be a numpy array"
-
     # Validation: p0 and p1 live in the same dimensional space
     assert len(p0) == len(p1), "Dimension mismatch. len(p0) != len(p1)"
 
@@ -47,7 +45,7 @@ def getIndecesBlock(p0,p1,ccMesh):
 
         ind  = np.where(indX & indY)
 
-    else:
+    elif dimMesh == 3:
         # Define the points
         x1 = p0[0]
         y1 = p0[1]
@@ -68,9 +66,9 @@ def getIndecesBlock(p0,p1,ccMesh):
 
 def defineBlockConductivity(p0,p1,ccMesh,condVals):
     """
-    Build a block with the conductivity specified by condVal.  Returns an array.
-    condVals[0]  conductivity of the block
-    condVals[1]  conductivity of the ground
+        Build a block with the conductivity specified by condVal.  Returns an array.
+        condVals[0]  conductivity of the block
+        condVals[1]  conductivity of the ground
     """
     sigma = np.zeros(ccMesh.shape[0]) + condVals[1]
     ind   = getIndecesBlock(p0,p1,ccMesh)
@@ -84,7 +82,8 @@ def defineTwoLayeredConductivity(depth,ccMesh,condVals):
     Define a two layered model.  Depth of the first layer must be specified.
     CondVals vector with the conductivity values of the layers.  Eg:
 
-    Convention to number the layers:
+    Convention to number the layers::
+
         <----------------------------|------------------------------------>
         0                          depth                                 zf
              1st layer                       2nd layer
@@ -98,13 +97,16 @@ def defineTwoLayeredConductivity(depth,ccMesh,condVals):
 
     # Identify 1st cell centered reference point
     p0[0] = ccMesh[0,0]
-    p0[1] = ccMesh[0,1]
-    p0[2] = ccMesh[0,2]
+    if dim>1: p0[1] = ccMesh[0,1]
+    if dim>2: p0[2] = ccMesh[0,2]
 
     # Identify the last cell-centered reference point
     p1[0] = ccMesh[-1,0]
-    p1[1] = ccMesh[-1,1]
-    p1[2] = ccMesh[-1,2] - depth;
+    if dim>1: p1[1] = ccMesh[-1,1]
+    if dim>2: p1[2] = ccMesh[-1,2]
+
+    # The depth is always defined on the last one.
+    p1[len(p1)-1] -= depth
 
     ind   = getIndecesBlock(p0,p1,ccMesh)
 
@@ -117,23 +119,24 @@ def scalarConductivity(ccMesh,pFunction):
     Define the distribution conductivity in the mesh according to the
     analytical expression given in pFunction
     """
-    xCC = ccMesh[:,0]
-    yCC = ccMesh[:,1]
-    zCC = ccMesh[:,2]
+    dim = np.size(ccMesh[0,:])
+    CC = [ccMesh[:,0]]
+    if dim>1: CC.append(ccMesh[:,1])
+    if dim>2: CC.append(ccMesh[:,2])
 
-    sigma = pFunction(xCC,yCC,zCC)
+
+    sigma = pFunction(*CC)
 
     return sigma
 
 if __name__ == '__main__':
 
-    import sys
-    sys.path.append('../')
-    from TensorMesh import TensorMesh
+    from SimPEG import TensorMesh
+    from matplotlib import pyplot as plt
 
     # Define the mesh
 
-    testDim = 3
+    testDim = 2
     h1 = 0.3*np.ones(7)
     h1[0] = 0.5
     h1[-1] = 0.6
@@ -157,8 +160,8 @@ if __name__ == '__main__':
     # ------------------- Test conductivities! --------------------------
     print('Testing 1 block conductivity')
 
-    p0 = np.array([0.5,0.5,0.5])
-    p1 = np.array([1.0,1.0,1.0])
+    p0 = np.array([0.5,0.5,0.5])[:testDim]
+    p1 = np.array([1.0,1.0,1.0])[:testDim]
     condVals = np.array([100,1e-6])
 
     sigma = defineBlockConductivity(p0,p1,ccMesh,condVals)
@@ -167,6 +170,7 @@ if __name__ == '__main__':
     print sigma.shape
     M.plotImage(sigma)
     print 'Done with block! :)'
+    plt.show()
 
     # -----------------------------------------
     print('Testing the two layered model')
@@ -178,11 +182,17 @@ if __name__ == '__main__':
     M.plotImage(sigma)
     print sigma
     print 'layer model!'
+    plt.show()
 
     # -----------------------------------------
     print('Testing scalar conductivity')
 
-    pFunction = lambda x,y,z: np.exp(x+y+z)
+    if testDim == 1:
+        pFunction = lambda x: np.exp(x)
+    elif testDim == 2:
+        pFunction = lambda x,y: np.exp(x+y)
+    elif testDim == 3:
+        pFunction = lambda x,y,z: np.exp(x+y+z)
 
     sigma = scalarConductivity(ccMesh,pFunction)
 
@@ -190,5 +200,6 @@ if __name__ == '__main__':
     M.plotImage(sigma)
     print sigma
     print 'Scalar conductivity defined!'
+    plt.show()
 
     # -----------------------------------------
diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py
index bb46a4e9..7c976ce3 100644
--- a/SimPEG/utils/__init__.py
+++ b/SimPEG/utils/__init__.py
@@ -1,4 +1,7 @@
+import matutils
+import sputils
+import lomutils
+import ModelBuilder
 from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
 from sputils import spzeros, kron3, speye, sdiag
 from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird
-import ModelBuilder
diff --git a/SimPEG/utils/lomutils.py b/SimPEG/utils/lomutils.py
index cdd6e0de..9020ce12 100644
--- a/SimPEG/utils/lomutils.py
+++ b/SimPEG/utils/lomutils.py
@@ -31,19 +31,18 @@ def volTetra(xyz, A, B, C, D):
     """
     Returns the volume for tetrahedras volume specified by the indexes A to D.
 
+    :param numpy.array xyz: X,Y,Z vertex vector
+    :param numpy.array A,B,C,D: vert index of the tetrahedra
+    :rtype: numpy.array
+    :return: V, volume of the tetrahedra
 
-    Input:
-       xyz      - X,Y,Z vertex vector
-       A,B,C,D  - vert index of the tetrahedra
+    Algorithm http://en.wikipedia.org/wiki/Tetrahedron#Volume
 
-    Output:
-       V        - volume
+    .. math::
 
-    Algorithm:     http://en.wikipedia.org/wiki/Tetrahedron#Volume
+       V = {1 \over 3} A  h
 
-       V = 1/3 A * h
-
-       V = 1/6 | ( a - d ) o ( ( b - d ) X ( c - d ) ) |
+       V = {1 \over 6} | ( a - d ) \cdot ( ( b - d )  ( c - d ) ) |
 
     """
 
@@ -69,7 +68,7 @@ def indexCube(nodes, gridSize, n=None):
     Output:
        index  - index in the order asked e.g. 'ABCD' --> (A,B,C,D)
 
-      TWO DIMENSIONS:
+    TWO DIMENSIONS::
 
       node(i,j)          node(i,j+1)
            A -------------- B
@@ -81,7 +80,7 @@ def indexCube(nodes, gridSize, n=None):
       node(i+1,j)        node(i+1,j+1)
 
 
-      THREE DIMENSIONS:
+    THREE DIMENSIONS::
 
             node(i,j,k+1)       node(i,j+1,k+1)
                 E --------------- F
@@ -97,10 +96,6 @@ def indexCube(nodes, gridSize, n=None):
            D -------------- C
       node(i+1,j,k)      node(i+1,j+1,k)
 
-
-         @author Rowan Cockett
-
-         Last modified on: 2013/07/26
     """
 
     assert type(nodes) == str, "Nodes must be a str variable: e.g. 'ABCD'"
@@ -211,7 +206,7 @@ def faceInfo(xyz, A, B, C, D, average=True, normalizeNormals=True):
     #
     # So also could be viewed as the average parallelogram.
     #
-    # WARNING: This does not compute correctly for concave quadrilaterals
+    # TODO: This does not compute correctly for concave quadrilaterals
     area = (length(nA)+length(nB)+length(nC)+length(nD))/4
 
     return N, area
diff --git a/SimPEG/utils/matutils.py b/SimPEG/utils/matutils.py
index 20cc33df..24286cdd 100644
--- a/SimPEG/utils/matutils.py
+++ b/SimPEG/utils/matutils.py
@@ -4,7 +4,7 @@ import numpy as np
 def mkvc(x, numDims=1):
     """Creates a vector with the number of dimension specified
 
-    e.g.:
+    e.g.::
 
         a = np.array([1, 2, 3])
 
@@ -43,7 +43,7 @@ def ndgrid(*args, **kwargs):
 
     The inputs can be a list or separate arguments.
 
-    e.g.
+    e.g.::
 
         a = np.array([1, 2, 3])
         b = np.array([1, 2])
diff --git a/docs/api_GaussNewton.rst b/docs/api_GaussNewton.rst
deleted file mode 100644
index 65012347..00000000
--- a/docs/api_GaussNewton.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-.. _api_GaussNewton:
-
-Gauss Newton
-************
-
-.. automodule:: SimPEG.GaussNewton
-    :members:
-    :undoc-members:
diff --git a/docs/api_Optimize.rst b/docs/api_Optimize.rst
new file mode 100644
index 00000000..9765bd3b
--- /dev/null
+++ b/docs/api_Optimize.rst
@@ -0,0 +1,8 @@
+.. _api_Inverse:
+
+Optimize
+********
+
+.. automodule:: SimPEG.inverse.Optimize
+    :members:
+    :undoc-members:
diff --git a/docs/api_Problem.rst b/docs/api_Problem.rst
new file mode 100644
index 00000000..83250f3e
--- /dev/null
+++ b/docs/api_Problem.rst
@@ -0,0 +1,26 @@
+.. _api_Problem:
+
+
+
+Problem
+*******
+
+.. automodule:: SimPEG.forward.Problem
+    :members:
+    :undoc-members:
+
+
+DCProblem
+*********
+
+.. automodule:: SimPEG.forward.DCProblem.DCProblem
+    :members:
+    :undoc-members:
+
+
+DCutils
+*******
+
+.. automodule:: SimPEG.forward.DCProblem.DCutils
+    :members:
+    :undoc-members:
diff --git a/docs/api_Tests.rst b/docs/api_Tests.rst
new file mode 100644
index 00000000..164f5b5d
--- /dev/null
+++ b/docs/api_Tests.rst
@@ -0,0 +1,8 @@
+.. _api_Tests:
+
+Testing SimPEG
+**************
+
+.. automodule:: SimPEG.tests.TestUtils
+    :members:
+    :undoc-members:
diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst
new file mode 100644
index 00000000..e952010b
--- /dev/null
+++ b/docs/api_Utils.rst
@@ -0,0 +1,21 @@
+.. _api_Utils:
+
+Utilities
+*********
+
+.. automodule:: SimPEG.utils.matutils
+    :members:
+    :undoc-members:
+
+.. automodule:: SimPEG.utils.sputils
+    :members:
+    :undoc-members:
+
+.. automodule:: SimPEG.utils.lomutils
+    :members:
+    :undoc-members:
+
+.. automodule:: SimPEG.utils.ModelBuilder
+    :members:
+    :undoc-members:
+
diff --git a/docs/index.rst b/docs/index.rst
index 22c9d5e9..4c01c4cd 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -30,20 +30,38 @@ Meshing & Operators
    api_DiffOperators
    api_InnerProducts
 
+Forward Problems
+================
+
+.. toctree::
+   :maxdepth: 2
+
+   api_Problem
+
 Inversion
 =========
 
 .. toctree::
    :maxdepth: 2
 
-   api_GaussNewton
+   api_Optimize
 
-Example Problems
-================
+Testing SimPEG
+==============
 
 .. toctree::
    :maxdepth: 2
 
+   api_Tests
+
+
+Utility Codes
+=============
+
+.. toctree::
+   :maxdepth: 2
+
+   api_Utils
 
 
 Project Index & Search