diff --git a/SimPEG/EM/TDEM/TDEM.py b/SimPEG/EM/TDEM/TDEM.py
index 7c451dc6..78f11617 100644
--- a/SimPEG/EM/TDEM/TDEM.py
+++ b/SimPEG/EM/TDEM/TDEM.py
@@ -62,7 +62,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
         return F
 
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         """
         Jvec computes the sensitivity times a vector
 
@@ -75,8 +75,8 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
             \mathbf{A} \\frac{d\mathbf{u}}{d\mathbf{m}} + \\frac{d\mathbf{A}(\mathbf{u})}{d\mathbf{m}} = \\frac{d \mathbf{RHS}}{d \mathbf{m}}
         """
 
-        if u is None:
-           u = self.fields(m)
+        if f is None:
+           f = self.fields(m)
 
         ftype = self._fieldType + 'Solution' # the thing we solved for
         self.curModel = m
@@ -115,17 +115,17 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
 
                 # here, we are lagging by a timestep, so filling in as we go
                 for projField in set([rx.projField for rx in src.rxList]):
-                    df_dmFun = getattr(u, '_%sDeriv'%projField, None)
+                    df_dmFun = getattr(f, '_%sDeriv'%projField, None)
                     # df_dm_v is dense, but we only need the times at (rx.P.T * ones > 0)
                     # This should be called rx.footprint
                     df_dm_v[src, '%sDeriv'%projField , tInd] = df_dmFun(tInd, src, dun_dm_v[:,i], v)
 
-                un_src = u[src,ftype,tInd+1]
+                un_src = f[src,ftype,tInd+1]
 
                 dA_dm_v   = self.getAdiagDeriv(tInd, un_src, v) # cell centered on time mesh
                 dRHS_dm_v = self.getRHSDeriv(tInd+1, src, v) # on nodes of time mesh
 
-                dAsubdiag_dm_v = self.getAsubdiagDeriv(tInd, u[src,ftype,tInd], v)
+                dAsubdiag_dm_v = self.getAsubdiagDeriv(tInd, f[src,ftype,tInd], v)
 
                 JRHS = dRHS_dm_v - dAsubdiag_dm_v - dA_dm_v
 
@@ -141,7 +141,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
         return Utils.mkvc(Jv)
 
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
 
         """
         Jvec computes the adjoint of the sensitivity times a vector
@@ -155,8 +155,8 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
             \\frac{d\mathbf{u}}{d\mathbf{m}} ^\\top \mathbf{A}^\\top  + \\frac{d\mathbf{A}(\mathbf{u})}{d\mathbf{m}} ^ \\top = \\frac{d \mathbf{RHS}}{d \mathbf{m}} ^ \\top
         """
 
-        if u is None:
-            u = self.fields(m)
+        if f is None:
+            f = self.fields(m)
 
         self.curModel = m
         ftype = self._fieldType + 'Solution' # the thing we solved for
@@ -166,7 +166,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
             v = self.dataPair(self.survey, v)
 
         df_duT_v = Fields_Derivs(self.mesh, self.survey)
-        ATinv_df_duT_v = np.zeros((len(self.survey.srcList), len(u[self.survey.srcList[0],ftype,0]))) # same size as fields at a single timestep
+        ATinv_df_duT_v = np.zeros((len(self.survey.srcList), len(f[self.survey.srcList[0],ftype,0]))) # same size as fields at a single timestep
 
         JTv = np.zeros(m.shape)
 
@@ -174,20 +174,20 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
         for src in self.survey.srcList:
             PT_v = Fields_Derivs(self.mesh, self.survey) # initialize storage for PT_v (don't need to preserve over sources)
             # initialize size
-            df_duT_v[src, '%sDeriv'%self._fieldType, :] = np.zeros_like(u[src, self._fieldType, :])
+            df_duT_v[src, '%sDeriv'%self._fieldType, :] = np.zeros_like(f[src, self._fieldType, :])
 
             for rx in src.rxList:
                 PT_v[src,'%sDeriv'%rx.projField,:] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True) # this is +=
 
                 # PT_v = np.reshape(curPT_v,(len(curPT_v)/self.timeMesh.nN, self.timeMesh.nN), order='F')
 
-                df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None)
+                df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
 
                 for tInd in range(self.nT+1):
                     cur = df_duTFun(tInd, src, None, Utils.mkvc(PT_v[src,'%sDeriv'%rx.projField,tInd]), adjoint=True)
 
                     df_duT_v[src, '%sDeriv'%self._fieldType, tInd] = df_duT_v[src, '%sDeriv'%self._fieldType, tInd] + Utils.mkvc(cur[0],2)
-                    JTv = JTv + cur[1]
+                    JTv = cur[1] + JTv
 
         del PT_v # no longer need this
 
@@ -227,10 +227,10 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
                     # (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
 
                 if tInd < self.nT - 1:
-                    dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, u[src,ftype,tInd+1], ATinv_df_duT_v[isrc,:], adjoint = True)
+                    dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, f[src,ftype,tInd+1], ATinv_df_duT_v[isrc,:], adjoint = True)
 
                 if tInd > -1:
-                    un_src = u[src,ftype,tInd+1]
+                    un_src = f[src,ftype,tInd+1]
                     dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh
                     dRHST_dm_v = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh
 
@@ -240,7 +240,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
                     # print np.linalg.norm(self.getInitialFieldsDeriv(src, df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1], adjoint=True))
                     # print np.linalg.norm(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1])
                     # vec = - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]) + Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1])
-                    # dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, u[src,ftype,tInd+1], Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint = True)
+                    # dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, f[src,ftype,tInd+1], Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint = True)
                     dRHST_dm_v = Utils.mkvc(self.getInitialFieldsDeriv(src, Utils.mkvc(ATinv_df_duT_v[isrc,:]) , adjoint=True))
 
                     JTv = JTv + Utils.mkvc( -dAsubdiagT_dm_v + dRHST_dm_v) #
@@ -253,7 +253,7 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
                     # JTv = JTv + Utils.mkvc(dRHST_dm_v0 + dRHST_dm_v1)
 
                     # print 'here'
-                    # inFields = self.getInitialFieldsDeriv(u[src,ftype,tInd+1], Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True)
+                    # inFields = self.getInitialFieldsDeriv(f[src,ftype,tInd+1], Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True)
                     #  # - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint=True)
                     # print inFields.shape
                     # JTv = JTv + inFields
@@ -268,11 +268,11 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
         # for src in self.survey.srcList:
         #     for projField in set(rx.projField):
         #         v = AdiagTinv * (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
-        #         JTv = JTv - Utils.mkvc(self.getAdiagDeriv(0, u[src, ftype, tInd], v, adjoint = True))
+        #         JTv = JTv - Utils.mkvc(self.getAdiagDeriv(0, f[src, ftype, tInd], v, adjoint = True))
         #         # JTv = JTv + self.getInitialFieldsDeriv(Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0] - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])), adjoint=True)
 
 
-        return Utils.mkvc(JTv)
+        return Utils.mkvc(JTv).astype(float)