diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py
index 7e7e548b..82b7b1e8 100644
--- a/conmech/solvers/optimization/optimization.py
+++ b/conmech/solvers/optimization/optimization.py
@@ -45,6 +45,7 @@
     "qsm",
     "qsmlm",
 }
+QSMLM_NAMES = {"dc " + name for name in QSMLM_NAMES}.union(QSMLM_NAMES)
 GLOBAL_QSMLM_NAMES = {
     "global quasi secant method",
     "global limited memory quasi secant method",
@@ -52,6 +53,7 @@
     "globqsm",
     "globqsmlm",
 }
+GLOBAL_QSMLM_NAMES = {"dc " + name for name in GLOBAL_QSMLM_NAMES}.union(GLOBAL_QSMLM_NAMES)
 
 
 class Optimization(Solver):
@@ -157,7 +159,11 @@ def _solve_impl(
             # pylint: disable=import-outside-toplevel,import-error)
             from kosopt.qsmlm import make_minimizer
 
-            self.minimizer = make_minimizer(self.loss, self.subgradient, self.sub2gradient)
+            self.minimizer = make_minimizer(
+                self.loss,
+                self.subgradient,
+                self.sub2gradient if method.lower().startswith("dc") else None
+            )
 
         while norm >= fixed_point_abs_tol:
             if method.lower() in QSMLM_NAMES:
diff --git a/examples/Makela_et_al_1998.py b/examples/Makela_et_al_1998.py
index dd3841d5..4f0806e3 100644
--- a/examples/Makela_et_al_1998.py
+++ b/examples/Makela_et_al_1998.py
@@ -247,11 +247,15 @@ def plot_losses(path):
     with open(path, "rb") as output:
         losses = pickle.load(output)
 
-    for mtd, values in losses.items():
+    for i, (mtd, values) in enumerate(losses.items()):
         forces_ = np.asarray(list(values.keys()))
         values_ = np.asarray(list(values.values()))[:, 0]
         times_ = np.asarray(list(values.values()))[:, 1]
-        plt.plot(forces_ / 1e3, -1 * values_, "-o", label=mtd)
+        # plt.plot(forces_ / 1e3, -1 * values_, "-o", label=mtd, linewidth=3*(len(losses.items()) - i))
+        total_width = 300
+        width = total_width / len(losses.items())
+        shift = total_width / 2 + i * width + width / 2
+        plt.bar(forces_ / 1e3 + shift, -1 * values_, label=mtd, width=width)
     plt.legend()
     plt.ylabel("$-\mathcal{L}(u)$")
     plt.xlabel(r"Load [kN/m$^2$]")
@@ -299,14 +303,14 @@ def loss_value(state, runner) -> float:
         Y[i] = MMLV99.subderivative_normal_direction(X[i], 0, 0)
     plt.plot(X, Y)
     plt.show()
-    # for i in range(1000):
-    #     Y[i] = MMLV99.sub2derivative_normal_direction(X[i], 0, 0)
-    # plt.plot(X, Y)
-    # plt.show()
-    # for i in range(1000):
-    #     Y[i] = MMLV99.subderivative_normal_direction(X[i], 0, 0) + MMLV99.sub2derivative_normal_direction(X[i], 0, 0)
-    # plt.plot(X, Y)
-    # plt.show()
+    for i in range(1000):
+        Y[i] = MMLV99.sub2derivative_normal_direction(X[i], 0, 0)
+    plt.plot(X, Y)
+    plt.show()
+    for i in range(1000):
+        Y[i] = MMLV99.subderivative_normal_direction(X[i], 0, 0) + MMLV99.sub2derivative_normal_direction(X[i], 0, 0)
+    plt.plot(X, Y)
+    plt.show()
 
     forces = np.asarray(
         (