diff --git a/moorpy/point.py b/moorpy/point.py
index 637ca10..358439a 100644
--- a/moorpy/point.py
+++ b/moorpy/point.py
@@ -341,7 +341,7 @@ def getStiffnessA(self, lines_only=False, xyz=False):
         else:                       # otherwise only return rows/columns of active DOFs
             return K[:,self.DOFs][self.DOFs,:]
 
-    def getDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, cbot, xyz=False):
+    def getDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, cbot):
         '''Gets inertia, added mass, damping, and stiffness matrices of Point due only to mooring lines (with other objects fixed)
         using a lumped mass model.
 
@@ -355,15 +355,19 @@ def getDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, cbot,
         #print("Getting Point "+str(self.number)+" analytic stiffness matrix...")
         
         def lump_matrix(matrix, endB):
-            if endB == 1:
-                matrix_woAnc = matrix[3:, 3:]
-                matrix_inv = np.linalg.inv(matrix_woAnc)
+            # The matrices should be symmetrical, but they can be slightly off due to numerical errors.
+            # Because we are going to invert them twice, we force them to be symmetrical to avoid amplifying the errors. 
+            matrix = (matrix + matrix.T)/2
+            if endB == 1:            
+                matrix_woAnc = matrix[3:, 3:] # For some reason that I don't know yet, we need to remove the elements in the opposite side of the matrix. Otherwise it doesn't work
+                matrix_inv = np.linalg.pinv(matrix_woAnc)
                 matrix_inv_coupled = matrix_inv[-3:, -3:]
             else:
                 matrix_woAnc = matrix[:-3, :-3]
-                matrix_inv = np.linalg.inv(matrix_woAnc)
+                matrix_inv = np.linalg.pinv(matrix_woAnc)
                 matrix_inv_coupled = matrix_inv[0:3, 0:3]
-            return np.linalg.inv(matrix_inv_coupled)
+            matrix_out = np.linalg.inv(matrix_inv_coupled)
+            return matrix_out 
 
         M = np.zeros([3,3])
         A = np.zeros([3,3])
@@ -382,11 +386,7 @@ def lump_matrix(matrix, endB):
             B += lump_matrix(B_all, endB)
             K += lump_matrix(K_all, endB)
 
-        if sum(np.isnan(K).ravel()) > 0: breakpoint()
-        if xyz:                     # if asked to output all DOFs, do it
-            return M, A, B, K
-        else:                       # otherwise only return rows/columns of active DOFs            
-            return M[:,self.DOFs][self.DOFs,:], A[:,self.DOFs][self.DOFs,:], B[:,self.DOFs][self.DOFs,:], K[:,self.DOFs][self.DOFs,:]
+        return M, A, B, K
         
     
     def getCost(self):
diff --git a/moorpy/system.py b/moorpy/system.py
index 90570f5..e831c02 100644
--- a/moorpy/system.py
+++ b/moorpy/system.py
@@ -3065,40 +3065,7 @@ def getSystemDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot,
                 A[i:i+n,i:i+n] += A1
                 B[i:i+n,i:i+n] += B1 
                 K[i:i+n,i:i+n] += K1
-                
-                # # go through attached lines and add cross-coupling terms
-                for lineID in point.attached:
-                    
-                    j = i + n
-                    
-                    # go through movable points to see if one is attached
-                    for point2 in self.pointList[self.pointList.index(point)+1: ]:
-                        if point2.type in d:
-                            if lineID in point2.attached:  # if this point is at the other end of the line
-
-                                # get cross-coupling stiffness of line: force on end attached to point1 due to motion of other end
-                                MB, AB, BB, KB = point2.getDynamicMatrices(omegas, S_zeta, r_dynamic_init, depth, kbot, cbot)
-                                if point.attachedEndB != 1:
-                                    MB = MB.T
-                                    AB = AB.T
-                                    BB = BB.T
-                                    KB = KB.T
-                                 
-                                # Trim stiffness matrix to only use the enabled DOFs of each point
-                                MB = MB[point.DOFs,:][:,point2.DOFs]
-                                AB = AB[point.DOFs,:][:,point2.DOFs]
-                                BB = BB[point.DOFs,:][:,point2.DOFs]
-                                KB = KB[point.DOFs,:][:,point2.DOFs]                                
-
-                                M[i:i+n          , j:j+point2.nDOF] += MB    # force on P1 due to movement of P2
-                                M[j:j+point2.nDOF, i:i+n          ] += MB.T  # mirror (f on P2 due to x of P1)
-                                A[i:i+n          , j:j+point2.nDOF] += AB  
-                                A[j:j+point2.nDOF, i:i+n          ] += AB.T
-                                B[i:i+n          , j:j+point2.nDOF] += BB  
-                                B[j:j+point2.nDOF, i:i+n          ] += BB.T
-                                K[i:i+n          , j:j+point2.nDOF] += KB  
-                                K[j:j+point2.nDOF, i:i+n          ] += KB.T                                
-                            j += point2.nDOF                  # if this point has DOFs we're considering, then count them                                
+                                             
                 i += n
         
         return M, A, B, K