--- /dev/null
+## The gradient should be zero at the optimal points.
+
+unit_test_equals("himmelblau_gradient(3,2) == 0", ...
+ 0, ...
+ himmelblau_gradient(3,2));
+
+x1 = -2.805118086952745;
+x2 = 3.131312518250573;
+msg = sprintf("himmelblau_gradient(%f, %f) == 0", x1, x2);
+unit_test_equals(msg, ...
+ true, ...
+ norm(himmelblau_gradient(x1, x2)) < 1e-13);
+
+x1 = -3.779310253377747;
+x2 = -3.283185991286170;
+msg = sprintf("himmelblau_gradient(%f, %f) == 0", x1, x2);
+unit_test_equals(msg, ...
+ true, ...
+ norm(himmelblau_gradient(x1, x2)) < 1e-13);
+
+x1 = 3.584428340330492;
+x2 = -1.848126526964404;
+msg = sprintf("himmelblau_gradient(%f, %f) == 0", x1, x2);
+unit_test_equals(msg, ...
+ true, ...
+ norm(himmelblau_gradient(x1, x2)) < 1e-13);