diff --git a/6_inv_free/NeoHookeanEnergy.py b/6_inv_free/NeoHookeanEnergy.py
index 9c13469..ec24dfa 100644
--- a/6_inv_free/NeoHookeanEnergy.py
+++ b/6_inv_free/NeoHookeanEnergy.py
@@ -72,10 +72,10 @@ def d2Psi_div_dF2(F, mu, lam):
                     rs = s * 2 + r
                     dP_div_dF[ij, rs] = M[0, 0] * U[i, 0] * VT[0, j] * U[r, 0] * VT[0, s] \
                         + M[0, 3] * U[i, 0] * VT[0, j] * U[r, 1] * VT[1, s] \
-                        + M[1, 1] * U[i, 0] * VT[1, j] * U[r, 0] * VT[1, s] \
-                        + M[1, 2] * U[i, 0] * VT[1, j] * U[r, 1] * VT[0, s] \
-                        + M[2, 1] * U[i, 1] * VT[0, j] * U[r, 0] * VT[1, s] \
-                        + M[2, 2] * U[i, 1] * VT[0, j] * U[r, 1] * VT[0, s] \
+                        + M[1, 1] * U[i, 1] * VT[0, j] * U[r, 1] * VT[0, s] \
+                        + M[1, 2] * U[i, 1] * VT[0, j] * U[r, 0] * VT[1, s] \
+                        + M[2, 1] * U[i, 0] * VT[1, j] * U[r, 1] * VT[0, s] \
+                        + M[2, 2] * U[i, 0] * VT[1, j] * U[r, 0] * VT[1, s] \
                         + M[3, 0] * U[i, 1] * VT[1, j] * U[r, 0] * VT[0, s] \
                         + M[3, 3] * U[i, 1] * VT[1, j] * U[r, 1] * VT[1, s]
     return dP_div_dF
diff --git a/7_self_contact/NeoHookeanEnergy.py b/7_self_contact/NeoHookeanEnergy.py
index d7a90a3..01d21e2 100644
--- a/7_self_contact/NeoHookeanEnergy.py
+++ b/7_self_contact/NeoHookeanEnergy.py
@@ -71,10 +71,10 @@ def d2Psi_div_dF2(F, mu, lam):
                     rs = s * 2 + r
                     dP_div_dF[ij, rs] = M[0, 0] * U[i, 0] * VT[0, j] * U[r, 0] * VT[0, s] \
                         + M[0, 3] * U[i, 0] * VT[0, j] * U[r, 1] * VT[1, s] \
-                        + M[1, 1] * U[i, 0] * VT[1, j] * U[r, 0] * VT[1, s] \
-                        + M[1, 2] * U[i, 0] * VT[1, j] * U[r, 1] * VT[0, s] \
-                        + M[2, 1] * U[i, 1] * VT[0, j] * U[r, 0] * VT[1, s] \
-                        + M[2, 2] * U[i, 1] * VT[0, j] * U[r, 1] * VT[0, s] \
+                        + M[1, 1] * U[i, 1] * VT[0, j] * U[r, 1] * VT[0, s] \
+                        + M[1, 2] * U[i, 1] * VT[0, j] * U[r, 0] * VT[1, s] \
+                        + M[2, 1] * U[i, 0] * VT[1, j] * U[r, 1] * VT[0, s] \
+                        + M[2, 2] * U[i, 0] * VT[1, j] * U[r, 0] * VT[1, s] \
                         + M[3, 0] * U[i, 1] * VT[1, j] * U[r, 0] * VT[0, s] \
                         + M[3, 3] * U[i, 1] * VT[1, j] * U[r, 1] * VT[1, s]
     return dP_div_dF
diff --git a/8_self_friction/NeoHookeanEnergy.py b/8_self_friction/NeoHookeanEnergy.py
index d7a90a3..01d21e2 100644
--- a/8_self_friction/NeoHookeanEnergy.py
+++ b/8_self_friction/NeoHookeanEnergy.py
@@ -71,10 +71,10 @@ def d2Psi_div_dF2(F, mu, lam):
                     rs = s * 2 + r
                     dP_div_dF[ij, rs] = M[0, 0] * U[i, 0] * VT[0, j] * U[r, 0] * VT[0, s] \
                         + M[0, 3] * U[i, 0] * VT[0, j] * U[r, 1] * VT[1, s] \
-                        + M[1, 1] * U[i, 0] * VT[1, j] * U[r, 0] * VT[1, s] \
-                        + M[1, 2] * U[i, 0] * VT[1, j] * U[r, 1] * VT[0, s] \
-                        + M[2, 1] * U[i, 1] * VT[0, j] * U[r, 0] * VT[1, s] \
-                        + M[2, 2] * U[i, 1] * VT[0, j] * U[r, 1] * VT[0, s] \
+                        + M[1, 1] * U[i, 1] * VT[0, j] * U[r, 1] * VT[0, s] \
+                        + M[1, 2] * U[i, 1] * VT[0, j] * U[r, 0] * VT[1, s] \
+                        + M[2, 1] * U[i, 0] * VT[1, j] * U[r, 1] * VT[0, s] \
+                        + M[2, 2] * U[i, 0] * VT[1, j] * U[r, 0] * VT[1, s] \
                         + M[3, 0] * U[i, 1] * VT[1, j] * U[r, 0] * VT[0, s] \
                         + M[3, 3] * U[i, 1] * VT[1, j] * U[r, 1] * VT[1, s]
     return dP_div_dF