diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py
index 16571fc3..e6eb38ce 100644
--- a/simpegEM/FDEM/FDEM.py
+++ b/simpegEM/FDEM/FDEM.py
@@ -28,8 +28,8 @@ class BaseFDEMProblem(BaseEMProblem):
         for freq in self.survey.freqs:
             A = self.getA(freq)
             rhs = RHS(freq)
-            solver = self.Solver(A, **self.solverOpts)
-            sol = solver.solve(rhs)
+            Ainv = self.Solver(A, **self.solverOpts)
+            sol = Ainv * rhs
             for fieldType in self.storeTheseFields:
                 Txs = self.survey.getTransmitters(freq)
                 F[Txs, fieldType] = CalcFields(sol, freq, fieldType)
@@ -46,12 +46,12 @@ class BaseFDEMProblem(BaseEMProblem):
 
         for freq in self.survey.freqs:
             A = self.getA(freq)
-            solver = self.Solver(A, **self.solverOpts)
+            Ainv = self.Solver(A, **self.solverOpts)
 
             for tx in self.survey.getTransmitters(freq):
                 u_tx = u[tx, self.solType]
                 w = self.getADeriv(freq, u_tx, v)
-                Ainvw = solver.solve(w)
+                Ainvw = Ainv * w
                 for rx in tx.rxList:
                     fAinvw = self.calcFields(Ainvw, freq, rx.projField)
                     P = lambda v: rx.projectFieldsDeriv(tx, self.mesh, u, v)
@@ -78,7 +78,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
         for freq in self.survey.freqs:
             AT = self.getA(freq).T
-            solver = self.Solver(AT, **self.solverOpts)
+            ATinv = self.Solver(AT, **self.solverOpts)
 
             for tx in self.survey.getTransmitters(freq):
                 u_tx = u[tx, self.solType]
@@ -87,7 +87,7 @@ class BaseFDEMProblem(BaseEMProblem):
                     PTv = rx.projectFieldsDeriv(tx, self.mesh, u, v[tx, rx], adjoint=True)
                     fPTv = self.calcFields(PTv, freq, rx.projField, adjoint=True)
 
-                    w = solver.solve( fPTv )
+                    w = ATinv * fPTv
                     Jtv_rx = - self.getADeriv(freq, u_tx, w, adjoint=True)
 
                     df_dm = self.calcFieldsDeriv(u_tx, freq, rx.projField, PTv, adjoint=True)
diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py
index 2782853d..c2d4e71b 100644
--- a/simpegEM/TDEM/BaseTDEM.py
+++ b/simpegEM/TDEM/BaseTDEM.py
@@ -32,22 +32,22 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         F = F or FieldsTDEM(self.mesh, self.survey)
 
         dtFact = None
-        Asolve = None
+        Ainv   = None
         for tInd, dt in enumerate(self.timeSteps):
             if dt != dtFact:
                 dtFact = dt
-                if Asolve is not None:
-                    Asolve.clean()
+                if Ainv is not None:
+                    Ainv.clean()
                 A = self.getA(tInd)
                 if self.verbose: print 'Factoring...   (dt = ' + str(dt) + ')'
-                Asolve = self.Solver(A, **self.solverOpts)
+                Ainv = self.Solver(A, **self.solverOpts)
                 if self.verbose: print 'Done'
             rhs = RHS(tInd, F)
-            sol = Asolve.solve(rhs)
+            sol = Ainv * rhs
             if sol.ndim == 1:
                 sol.shape = (sol.size,1)
             F[:,:,tInd+1] = CalcFields(sol, tInd)
-        Asolve.clean()
+        Ainv.clean()
         return F
 
     def adjoint(self, m, RHS, CalcFields, F=None):
@@ -55,22 +55,22 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         F = F or FieldsTDEM(self.mesh, self.survey)
 
         dtFact = None
-        Asolve = None
+        Ainv   = None
         for tInd, dt in reversed(list(enumerate(self.timeSteps))):
             if dt != dtFact:
                 dtFact = dt
-                if Asolve is not None:
-                    Asolve.clean()
+                if Ainv is not None:
+                    Ainv.clean()
                 A = self.getA(tInd)
                 if self.verbose: print 'Factoring...   (dt = ' + str(dt) + ')'
-                Asolve = self.Solver(A, **self.solverOpts)
+                Ainv = self.Solver(A, **self.solverOpts)
                 if self.verbose: print 'Done'
             rhs = RHS(tInd, F)
-            sol = Asolve.solve(rhs)
+            sol = Ainv * rhs
             if sol.ndim == 1:
                 sol.shape = (sol.size,1)
             F[:,:,tInd+1] = CalcFields(sol, tInd)
-        Asolve.clean()
+        Ainv.clean()
         return F
 
     def Jvec(self, m, v, u=None):