diff --git a/lppls/lppls.py b/lppls/lppls.py
index 54e7c63..186c471 100644
--- a/lppls/lppls.py
+++ b/lppls/lppls.py
@@ -27,7 +27,8 @@ def __init__(self, observations):
     @staticmethod
     @njit
     def lppls(t, tc, m, w, a, b, c1, c2):
-        return a + np.power(tc - t, m) * (b + ((c1 * np.cos(w * np.log(tc - t))) + (c2 * np.sin(w * np.log(tc - t)))))
+        dt = np.abs(tc - t) + 1e-8
+        return a + np.power(dt, m) * (b + ((c1 * np.cos(w * np.log(dt))) + (c2 * np.sin(w * np.log(dt)))))
 
     def func_restricted(self, x, *args):
         """
@@ -49,7 +50,7 @@ def func_restricted(self, x, *args):
         # print('type', type(res))
         # print('func_restricted', res)
 
-        delta = [self.lppls(t, tc, m, w, a, b, c1, c2) for t in observations[0, :]]
+        delta = self.lppls(observations[0, :], tc, m, w, a, b, c1, c2)
         delta = np.subtract(delta, observations[1, :])
         delta = np.power(delta, 2)
         return np.sum(delta)
@@ -64,8 +65,7 @@ def matrix_equation(observations, tc, m, w):
         P = observations[1]
         N = len(T)
 
-        # @TODO make taking tc - t or |tc - t| configurable
-        dT = np.abs(tc - T)
+        dT = np.abs(tc - T) + 1e-8
         phase = np.log(dT)
 
         fi = np.power(dT, m)
@@ -99,6 +99,8 @@ def matrix_equation(observations, tc, m, w):
             [np.sum(yihi)]
         ])
 
+        matrix_1 += 1e-8 * np.eye(matrix_1.shape[0])
+
         return np.linalg.solve(matrix_1, matrix_2)
 
     def fit(self, max_searches, minimizer='Nelder-Mead', obs=None):
diff --git a/setup.py b/setup.py
index 039e4fd..9706ee0 100644
--- a/setup.py
+++ b/setup.py
@@ -4,7 +4,7 @@
     long_description = fh.read()
 
 setuptools.setup(name='lppls',
-      version='0.6.12',
+      version='0.6.13',
       description='A Python module for fitting the LPPLS model to data.',
       packages=['lppls'],
       author='Josh Nielsen',