diff --git a/docs/source/gradient_tutorial.rst b/docs/source/gradient_tutorial.rst
index 3fef6b6..aea407a 100644
--- a/docs/source/gradient_tutorial.rst
+++ b/docs/source/gradient_tutorial.rst
@@ -8,46 +8,46 @@
 General Unconstrained Minimization
 ==================================
 
-While much of Ceres Solver is devoted to solving non-linear least
-squares problems, internally it contains a solver that can solve
-general unconstrained optimization problems using just their objective
-function value and gradients. The ``GradientProblem`` and
-``GradientProblemSolver`` objects give the user access to this solver.
-
-So without much further ado, let us look at how one goes about using
-them.
+Ceres Solver besides being able to solve non-linear least squares
+problem can also solve general unconstrained problems using just their
+objective function value and gradients. In this chapter we will see
+how to do this.
 
 Rosenbrock's Function
 =====================
 
-We consider the minimization of the famous `Rosenbrock's function
+Consider minimizing the famous `Rosenbrock's function
 <http://en.wikipedia.org/wiki/Rosenbrock_function>`_ [#f1]_.
 
-We begin by defining an instance of the ``FirstOrderFunction``
-interface. This is the object that is responsible for computing the
-objective function value and the gradient (if required). This is the
-analog of the :class:`CostFunction` when defining non-linear least
-squares problems in Ceres.
+The simplest way to minimize is to define a templated functor to
+evaluate the objective value of this function and then use Ceres
+Solver's automatic differentiation to compute its derivatives.
+
+We begin by defining a templated functor and then using
+``AutoDiffFirstOrderFunction`` to construct an instance of the
+``FirstOrderFunction`` interface. This is the object that is
+responsible for computing the objective function value and the
+gradient (if required). This is the analog of the
+:class:`CostFunction` when defining non-linear least squares problems
+in Ceres.
 
 .. code::
 
-  class Rosenbrock : public ceres::FirstOrderFunction {
-   public:
-    virtual bool Evaluate(const double* parameters,
-                          double* cost,
-                          double* gradient) const {
-      const double x = parameters[0];
-      const double y = parameters[1];
-
+  // f(x,y) = (1-x)^2 + 100(y - x^2)^2;
+  struct Rosenbrock {
+    template <typename T>
+    bool operator()(const T* parameters, T* cost) const {
+      const T x = parameters[0];
+      const T y = parameters[1];
       cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x);
-      if (gradient != nullptr) {
-        gradient[0] = -2.0 * (1.0 - x) - 200.0 * (y - x * x) * 2.0 * x;
-        gradient[1] = 200.0 * (y - x * x);
-      }
       return true;
     }
 
-    virtual int NumParameters() const { return 2; }
+    static ceres::FirstOrderFunction* Create() {
+      constexpr int kNumParameters = 2;
+      return new ceres::AutoDiffFirstOrderFunction<Rosenbrock, kNumParameters>(
+          new Rosenbrock);
+    }
   };
 
 
@@ -58,7 +58,7 @@
 
     double parameters[2] = {-1.2, 1.0};
 
-    ceres::GradientProblem problem(new Rosenbrock());
+    ceres::GradientProblem problem(Rosenbrock::Create());
 
     ceres::GradientProblemSolver::Options options;
     options.minimizer_progress_to_stdout = true;
@@ -74,43 +74,43 @@
 
 .. code-block:: bash
 
-     0: f: 2.420000e+01 d: 0.00e+00 g: 2.16e+02 h: 0.00e+00 s: 0.00e+00 e:  0 it: 2.00e-05 tt: 2.00e-05
-     1: f: 4.280493e+00 d: 1.99e+01 g: 1.52e+01 h: 2.01e-01 s: 8.62e-04 e:  2 it: 7.32e-05 tt: 2.19e-04
-     2: f: 3.571154e+00 d: 7.09e-01 g: 1.35e+01 h: 3.78e-01 s: 1.34e-01 e:  3 it: 2.50e-05 tt: 2.68e-04
-     3: f: 3.440869e+00 d: 1.30e-01 g: 1.73e+01 h: 1.36e-01 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 2.92e-04
-     4: f: 3.213597e+00 d: 2.27e-01 g: 1.55e+01 h: 1.06e-01 s: 4.59e-01 e:  1 it: 2.86e-06 tt: 3.14e-04
-     5: f: 2.839723e+00 d: 3.74e-01 g: 1.05e+01 h: 1.34e-01 s: 5.24e-01 e:  1 it: 2.86e-06 tt: 3.36e-04
-     6: f: 2.448490e+00 d: 3.91e-01 g: 1.29e+01 h: 3.04e-01 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 3.58e-04
-     7: f: 1.943019e+00 d: 5.05e-01 g: 4.00e+00 h: 8.81e-02 s: 7.43e-01 e:  1 it: 4.05e-06 tt: 3.79e-04
-     8: f: 1.731469e+00 d: 2.12e-01 g: 7.36e+00 h: 1.71e-01 s: 4.60e-01 e:  2 it: 9.06e-06 tt: 4.06e-04
-     9: f: 1.503267e+00 d: 2.28e-01 g: 6.47e+00 h: 8.66e-02 s: 1.00e+00 e:  1 it: 3.81e-06 tt: 4.33e-04
-    10: f: 1.228331e+00 d: 2.75e-01 g: 2.00e+00 h: 7.70e-02 s: 7.90e-01 e:  1 it: 3.81e-06 tt: 4.54e-04
-    11: f: 1.016523e+00 d: 2.12e-01 g: 5.15e+00 h: 1.39e-01 s: 3.76e-01 e:  2 it: 1.00e-05 tt: 4.82e-04
-    12: f: 9.145773e-01 d: 1.02e-01 g: 6.74e+00 h: 7.98e-02 s: 1.00e+00 e:  1 it: 3.10e-06 tt: 5.03e-04
-    13: f: 7.508302e-01 d: 1.64e-01 g: 3.88e+00 h: 5.76e-02 s: 4.93e-01 e:  1 it: 2.86e-06 tt: 5.25e-04
-    14: f: 5.832378e-01 d: 1.68e-01 g: 5.56e+00 h: 1.42e-01 s: 1.00e+00 e:  1 it: 3.81e-06 tt: 5.47e-04
-    15: f: 3.969581e-01 d: 1.86e-01 g: 1.64e+00 h: 1.17e-01 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 5.68e-04
-    16: f: 3.171557e-01 d: 7.98e-02 g: 3.84e+00 h: 1.18e-01 s: 3.97e-01 e:  2 it: 9.06e-06 tt: 5.94e-04
-    17: f: 2.641257e-01 d: 5.30e-02 g: 3.27e+00 h: 6.14e-02 s: 1.00e+00 e:  1 it: 3.10e-06 tt: 6.16e-04
-    18: f: 1.909730e-01 d: 7.32e-02 g: 5.29e-01 h: 8.55e-02 s: 6.82e-01 e:  1 it: 4.05e-06 tt: 6.42e-04
-    19: f: 1.472012e-01 d: 4.38e-02 g: 3.11e+00 h: 1.20e-01 s: 3.47e-01 e:  2 it: 1.00e-05 tt: 6.69e-04
-    20: f: 1.093558e-01 d: 3.78e-02 g: 2.97e+00 h: 8.43e-02 s: 1.00e+00 e:  1 it: 3.81e-06 tt: 6.91e-04
-    21: f: 6.710346e-02 d: 4.23e-02 g: 1.42e+00 h: 9.64e-02 s: 8.85e-01 e:  1 it: 3.81e-06 tt: 7.12e-04
-    22: f: 3.993377e-02 d: 2.72e-02 g: 2.30e+00 h: 1.29e-01 s: 4.63e-01 e:  2 it: 9.06e-06 tt: 7.39e-04
-    23: f: 2.911794e-02 d: 1.08e-02 g: 2.55e+00 h: 6.55e-02 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 7.62e-04
-    24: f: 1.457683e-02 d: 1.45e-02 g: 2.77e-01 h: 6.37e-02 s: 6.14e-01 e:  1 it: 3.81e-06 tt: 7.84e-04
-    25: f: 8.577515e-03 d: 6.00e-03 g: 2.86e+00 h: 1.40e-01 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 8.05e-04
-    26: f: 3.486574e-03 d: 5.09e-03 g: 1.76e-01 h: 1.23e-02 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 8.27e-04
-    27: f: 1.257570e-03 d: 2.23e-03 g: 1.39e-01 h: 5.08e-02 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 8.48e-04
-    28: f: 2.783568e-04 d: 9.79e-04 g: 6.20e-01 h: 6.47e-02 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 8.69e-04
-    29: f: 2.533399e-05 d: 2.53e-04 g: 1.68e-02 h: 1.98e-03 s: 1.00e+00 e:  1 it: 3.81e-06 tt: 8.91e-04
-    30: f: 7.591572e-07 d: 2.46e-05 g: 5.40e-03 h: 9.27e-03 s: 1.00e+00 e:  1 it: 3.81e-06 tt: 9.12e-04
-    31: f: 1.902460e-09 d: 7.57e-07 g: 1.62e-03 h: 1.89e-03 s: 1.00e+00 e:  1 it: 2.86e-06 tt: 9.33e-04
-    32: f: 1.003030e-12 d: 1.90e-09 g: 3.50e-05 h: 3.52e-05 s: 1.00e+00 e:  1 it: 3.10e-06 tt: 9.54e-04
-    33: f: 4.835994e-17 d: 1.00e-12 g: 1.05e-07 h: 1.13e-06 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 9.81e-04
-    34: f: 1.885250e-22 d: 4.84e-17 g: 2.69e-10 h: 1.45e-08 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 1.00e-03
+     0: f: 2.420000e+01 d: 0.00e+00 g: 2.16e+02 h: 0.00e+00 s: 0.00e+00 e:  0 it: 2.29e-05 tt: 2.29e-05
+     1: f: 4.280493e+00 d: 1.99e+01 g: 1.52e+01 h: 2.01e-01 s: 8.62e-04 e:  2 it: 8.39e-05 tt: 1.62e-04
+     2: f: 3.571154e+00 d: 7.09e-01 g: 1.35e+01 h: 3.78e-01 s: 1.34e-01 e:  3 it: 2.22e-05 tt: 1.91e-04
+     3: f: 3.440869e+00 d: 1.30e-01 g: 1.73e+01 h: 1.36e-01 s: 1.00e+00 e:  1 it: 5.01e-06 tt: 2.01e-04
+     4: f: 3.213597e+00 d: 2.27e-01 g: 1.55e+01 h: 1.06e-01 s: 4.59e-01 e:  1 it: 3.81e-06 tt: 2.10e-04
+     5: f: 2.839723e+00 d: 3.74e-01 g: 1.05e+01 h: 1.34e-01 s: 5.24e-01 e:  1 it: 4.05e-06 tt: 2.19e-04
+     6: f: 2.448490e+00 d: 3.91e-01 g: 1.29e+01 h: 3.04e-01 s: 1.00e+00 e:  1 it: 5.01e-06 tt: 2.28e-04
+     7: f: 1.943019e+00 d: 5.05e-01 g: 4.00e+00 h: 8.81e-02 s: 7.43e-01 e:  1 it: 4.05e-06 tt: 2.36e-04
+     8: f: 1.731469e+00 d: 2.12e-01 g: 7.36e+00 h: 1.71e-01 s: 4.60e-01 e:  2 it: 1.22e-05 tt: 2.52e-04
+     9: f: 1.503267e+00 d: 2.28e-01 g: 6.47e+00 h: 8.66e-02 s: 1.00e+00 e:  1 it: 5.96e-06 tt: 2.66e-04
+    10: f: 1.228331e+00 d: 2.75e-01 g: 2.00e+00 h: 7.70e-02 s: 7.90e-01 e:  1 it: 4.05e-06 tt: 2.75e-04
+    11: f: 1.016523e+00 d: 2.12e-01 g: 5.15e+00 h: 1.39e-01 s: 3.76e-01 e:  2 it: 9.06e-06 tt: 2.88e-04
+    12: f: 9.145773e-01 d: 1.02e-01 g: 6.74e+00 h: 7.98e-02 s: 1.00e+00 e:  1 it: 5.01e-06 tt: 2.97e-04
+    13: f: 7.508302e-01 d: 1.64e-01 g: 3.88e+00 h: 5.76e-02 s: 4.93e-01 e:  1 it: 5.01e-06 tt: 3.05e-04
+    14: f: 5.832378e-01 d: 1.68e-01 g: 5.56e+00 h: 1.42e-01 s: 1.00e+00 e:  1 it: 4.77e-06 tt: 3.13e-04
+    15: f: 3.969581e-01 d: 1.86e-01 g: 1.64e+00 h: 1.17e-01 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 3.20e-04
+    16: f: 3.171557e-01 d: 7.98e-02 g: 3.84e+00 h: 1.18e-01 s: 3.97e-01 e:  2 it: 8.82e-06 tt: 3.33e-04
+    17: f: 2.641257e-01 d: 5.30e-02 g: 3.27e+00 h: 6.14e-02 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 3.42e-04
+    18: f: 1.909730e-01 d: 7.32e-02 g: 5.29e-01 h: 8.55e-02 s: 6.82e-01 e:  1 it: 1.00e-05 tt: 4.64e-04
+    19: f: 1.472012e-01 d: 4.38e-02 g: 3.11e+00 h: 1.20e-01 s: 3.47e-01 e:  2 it: 1.29e-05 tt: 4.87e-04
+    20: f: 1.093558e-01 d: 3.78e-02 g: 2.97e+00 h: 8.43e-02 s: 1.00e+00 e:  1 it: 5.01e-06 tt: 4.97e-04
+    21: f: 6.710346e-02 d: 4.23e-02 g: 1.42e+00 h: 9.64e-02 s: 8.85e-01 e:  1 it: 4.05e-06 tt: 5.06e-04
+    22: f: 3.993377e-02 d: 2.72e-02 g: 2.30e+00 h: 1.29e-01 s: 4.63e-01 e:  2 it: 1.00e-05 tt: 5.25e-04
+    23: f: 2.911794e-02 d: 1.08e-02 g: 2.55e+00 h: 6.55e-02 s: 1.00e+00 e:  1 it: 5.01e-06 tt: 5.34e-04
+    24: f: 1.457683e-02 d: 1.45e-02 g: 2.77e-01 h: 6.37e-02 s: 6.14e-01 e:  1 it: 4.05e-06 tt: 5.42e-04
+    25: f: 8.577515e-03 d: 6.00e-03 g: 2.86e+00 h: 1.40e-01 s: 1.00e+00 e:  1 it: 3.81e-06 tt: 5.49e-04
+    26: f: 3.486574e-03 d: 5.09e-03 g: 1.76e-01 h: 1.23e-02 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 5.57e-04
+    27: f: 1.257570e-03 d: 2.23e-03 g: 1.39e-01 h: 5.08e-02 s: 1.00e+00 e:  1 it: 3.81e-06 tt: 5.65e-04
+    28: f: 2.783568e-04 d: 9.79e-04 g: 6.20e-01 h: 6.47e-02 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 5.73e-04
+    29: f: 2.533399e-05 d: 2.53e-04 g: 1.68e-02 h: 1.98e-03 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 5.81e-04
+    30: f: 7.591572e-07 d: 2.46e-05 g: 5.40e-03 h: 9.27e-03 s: 1.00e+00 e:  1 it: 5.96e-06 tt: 6.30e-04
+    31: f: 1.902460e-09 d: 7.57e-07 g: 1.62e-03 h: 1.89e-03 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 6.39e-04
+    32: f: 1.003030e-12 d: 1.90e-09 g: 3.50e-05 h: 3.52e-05 s: 1.00e+00 e:  1 it: 3.81e-06 tt: 6.47e-04
+    33: f: 4.835994e-17 d: 1.00e-12 g: 1.05e-07 h: 1.13e-06 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 6.59e-04
+    34: f: 1.885250e-22 d: 4.84e-17 g: 2.69e-10 h: 1.45e-08 s: 1.00e+00 e:  1 it: 4.05e-06 tt: 6.67e-04
 
-  Solver Summary (v 1.12.0-lapack-suitesparse-cxsparse-no_openmp)
+  Solver Summary (v 2.0.0-eigen-(3.3.9)-lapack-suitesparse-(5.8.1)-cxsparse-(3.2.0)-acceleratesparse-eigensparse-no_openmp-no_custom_blas)
 
   Parameters                                  2
   Line search direction              LBFGS (20)
@@ -119,20 +119,86 @@
 
   Cost:
   Initial                          2.420000e+01
-  Final                            1.885250e-22
+  Final                            1.955192e-27
   Change                           2.420000e+01
 
-  Minimizer iterations                       35
+  Minimizer iterations                       36
 
   Time (in seconds):
 
-    Cost evaluation                       0.000
-    Gradient evaluation                   0.000
-  Total                                   0.003
+    Cost evaluation                    0.000000 (0)
+    Gradient & cost evaluation         0.000008 (44)
+    Polynomial minimization            0.000067
+  Total                                0.000721
 
-  Termination:                      CONVERGENCE (Gradient tolerance reached. Gradient max norm: 9.032775e-13 <= 1.000000e-10)
+  Termination:                      CONVERGENCE (Parameter tolerance reached. Relative step_norm: 1.890726e-11 <= 1.000000e-08.)
+
+  Initial x: -1.2 y: 1
+  Final   x: 1 y: 1
+
+
+If you are unable to use automatic differentiation for some reason
+(say beacause you need to call an external library), then you can
+use numeric differentiation. In that case the functor is defined as
+follows [#f2]_.
+
+.. code::
+
+  // f(x,y) = (1-x)^2 + 100(y - x^2)^2;
+  struct Rosenbrock {
+    bool operator()(const double* parameters, double* cost) const {
+      const double x = parameters[0];
+      const double y = parameters[1];
+      cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x);
+      return true;
+    }
+
+    static ceres::FirstOrderFunction* Create() {
+      constexpr int kNumParameters = 2;
+      return new ceres::NumericDiffFirstOrderFunction<Rosenbrock,
+                                                      ceres::CENTRAL,
+                                                      kNumParameters>(
+          new Rosenbrock);
+    }
+};
+
+And finally, if you would rather compute the derivatives by hand (say
+because the size of the parameter vector is too large to be
+automatically differentiated). Then you should define an instance of
+``FirstOrderFunction``, which is the analog of :class:`CostFunction`
+for non-linear least squares problems [#f3]_.
+
+.. code::
+
+  // f(x,y) = (1-x)^2 + 100(y - x^2)^2;
+  class Rosenbrock final  : public ceres::FirstOrderFunction {
+    public:
+      ~Rosenbrock() override {}
+
+      bool Evaluate(const double* parameters,
+                             double* cost,
+                             double* gradient) const override {
+         const double x = parameters[0];
+         const double y = parameters[1];
+
+         cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x);
+         if (gradient) {
+           gradient[0] = -2.0 * (1.0 - x) - 200.0 * (y - x * x) * 2.0 * x;
+           gradient[1] = 200.0 * (y - x * x);
+         }
+        return true;
+     }
+
+     int NumParameters() const override { return 2; }
+  };
 
 .. rubric:: Footnotes
 
 .. [#f1] `examples/rosenbrock.cc
    <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/rosenbrock.cc>`_
+
+.. [#f2] `examples/rosenbrock_numeric_diff.cc
+   <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/rosenbrock_numeric_diff.cc>`_
+
+.. [#f3] `examples/rosenbrock_analytic_diff.cc
+   <https://ceres-solver.googlesource.com/ceres-solver/+/master/examples/rosenbrock_analytic_diff.cc>`_
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 7f9b117..a43c8aa 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -50,6 +50,12 @@
 add_executable(rosenbrock rosenbrock.cc)
 target_link_libraries(rosenbrock Ceres::ceres)
 
+add_executable(rosenbrock_analytic_diff rosenbrock_analytic_diff.cc)
+target_link_libraries(rosenbrock_analytic_diff Ceres::ceres)
+
+add_executable(rosenbrock_numeric_diff rosenbrock_numeric_diff.cc)
+target_link_libraries(rosenbrock_numeric_diff Ceres::ceres)
+
 add_executable(curve_fitting_c curve_fitting.c)
 target_link_libraries(curve_fitting_c Ceres::ceres)
 # Force CMake to link curve_fitting_c using the C linker.
diff --git a/examples/rosenbrock.cc b/examples/rosenbrock.cc
index 1b9aef6..8243457 100644
--- a/examples/rosenbrock.cc
+++ b/examples/rosenbrock.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2021 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -32,25 +32,20 @@
 #include "glog/logging.h"
 
 // f(x,y) = (1-x)^2 + 100(y - x^2)^2;
-class Rosenbrock : public ceres::FirstOrderFunction {
- public:
-  virtual ~Rosenbrock() {}
-
-  virtual bool Evaluate(const double* parameters,
-                        double* cost,
-                        double* gradient) const {
-    const double x = parameters[0];
-    const double y = parameters[1];
-
+struct Rosenbrock {
+  template <typename T>
+  bool operator()(const T* parameters, T* cost) const {
+    const T x = parameters[0];
+    const T y = parameters[1];
     cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x);
-    if (gradient != NULL) {
-      gradient[0] = -2.0 * (1.0 - x) - 200.0 * (y - x * x) * 2.0 * x;
-      gradient[1] = 200.0 * (y - x * x);
-    }
     return true;
   }
 
-  virtual int NumParameters() const { return 2; }
+  static ceres::FirstOrderFunction* Create() {
+    constexpr int kNumParameters = 2;
+    return new ceres::AutoDiffFirstOrderFunction<Rosenbrock, kNumParameters>(
+        new Rosenbrock);
+  }
 };
 
 int main(int argc, char** argv) {
@@ -62,7 +57,7 @@
   options.minimizer_progress_to_stdout = true;
 
   ceres::GradientProblemSolver::Summary summary;
-  ceres::GradientProblem problem(new Rosenbrock());
+  ceres::GradientProblem problem(Rosenbrock::Create());
   ceres::Solve(options, problem, parameters, &summary);
 
   std::cout << summary.FullReport() << "\n";
diff --git a/examples/rosenbrock_analytic_diff.cc b/examples/rosenbrock_analytic_diff.cc
new file mode 100644
index 0000000..b133f5a
--- /dev/null
+++ b/examples/rosenbrock_analytic_diff.cc
@@ -0,0 +1,75 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2021 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/ceres.h"
+#include "glog/logging.h"
+
+// f(x,y) = (1-x)^2 + 100(y - x^2)^2;
+class Rosenbrock final : public ceres::FirstOrderFunction {
+ public:
+  ~Rosenbrock() override {}
+
+  bool Evaluate(const double* parameters,
+                double* cost,
+                double* gradient) const override {
+    const double x = parameters[0];
+    const double y = parameters[1];
+
+    cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x);
+
+    if (gradient) {
+      gradient[0] = -2.0 * (1.0 - x) - 200.0 * (y - x * x) * 2.0 * x;
+      gradient[1] = 200.0 * (y - x * x);
+    }
+
+    return true;
+  }
+
+  int NumParameters() const override { return 2; }
+};
+
+int main(int argc, char** argv) {
+  google::InitGoogleLogging(argv[0]);
+
+  double parameters[2] = {-1.2, 1.0};
+
+  ceres::GradientProblemSolver::Options options;
+  options.minimizer_progress_to_stdout = true;
+
+  ceres::GradientProblemSolver::Summary summary;
+  ceres::GradientProblem problem(new Rosenbrock());
+  ceres::Solve(options, problem, parameters, &summary);
+
+  std::cout << summary.FullReport() << "\n";
+  std::cout << "Initial x: " << -1.2 << " y: " << 1.0 << "\n";
+  std::cout << "Final   x: " << parameters[0] << " y: " << parameters[1]
+            << "\n";
+  return 0;
+}
diff --git a/examples/rosenbrock_numeric_diff.cc b/examples/rosenbrock_numeric_diff.cc
new file mode 100644
index 0000000..f2f6e8a
--- /dev/null
+++ b/examples/rosenbrock_numeric_diff.cc
@@ -0,0 +1,69 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2021 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/ceres.h"
+#include "glog/logging.h"
+
+// f(x,y) = (1-x)^2 + 100(y - x^2)^2;
+struct Rosenbrock {
+  bool operator()(const double* parameters, double* cost) const {
+    const double x = parameters[0];
+    const double y = parameters[1];
+    cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x);
+    return true;
+  }
+
+  static ceres::FirstOrderFunction* Create() {
+    constexpr int kNumParameters = 2;
+    return new ceres::NumericDiffFirstOrderFunction<Rosenbrock,
+                                                    ceres::CENTRAL,
+                                                    kNumParameters>(
+        new Rosenbrock);
+  }
+};
+
+int main(int argc, char** argv) {
+  google::InitGoogleLogging(argv[0]);
+
+  double parameters[2] = {-1.2, 1.0};
+
+  ceres::GradientProblemSolver::Options options;
+  options.minimizer_progress_to_stdout = true;
+
+  ceres::GradientProblemSolver::Summary summary;
+  ceres::GradientProblem problem(Rosenbrock::Create());
+  ceres::Solve(options, problem, parameters, &summary);
+
+  std::cout << summary.FullReport() << "\n";
+  std::cout << "Initial x: " << -1.2 << " y: " << 1.0 << "\n";
+  std::cout << "Final   x: " << parameters[0] << " y: " << parameters[1]
+            << "\n";
+  return 0;
+}
diff --git a/include/ceres/ceres.h b/include/ceres/ceres.h
index d249351..86a84ab 100644
--- a/include/ceres/ceres.h
+++ b/include/ceres/ceres.h
@@ -35,6 +35,7 @@
 #define CERES_PUBLIC_CERES_H_
 
 #include "ceres/autodiff_cost_function.h"
+#include "ceres/autodiff_first_order_function.h"
 #include "ceres/autodiff_local_parameterization.h"
 #include "ceres/conditioned_cost_function.h"
 #include "ceres/context.h"
@@ -47,6 +48,7 @@
 #include "ceres/dynamic_cost_function_to_functor.h"
 #include "ceres/dynamic_numeric_diff_cost_function.h"
 #include "ceres/evaluation_callback.h"
+#include "ceres/first_order_function.h"
 #include "ceres/gradient_checker.h"
 #include "ceres/gradient_problem.h"
 #include "ceres/gradient_problem_solver.h"
@@ -55,6 +57,7 @@
 #include "ceres/local_parameterization.h"
 #include "ceres/loss_function.h"
 #include "ceres/numeric_diff_cost_function.h"
+#include "ceres/numeric_diff_first_order_function.h"
 #include "ceres/numeric_diff_options.h"
 #include "ceres/ordered_groups.h"
 #include "ceres/problem.h"
diff --git a/include/ceres/numeric_diff_first_order_function.h b/include/ceres/numeric_diff_first_order_function.h
new file mode 100644
index 0000000..8663835
--- /dev/null
+++ b/include/ceres/numeric_diff_first_order_function.h
@@ -0,0 +1,162 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2019 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#ifndef CERES_PUBLIC_NUMERIC_DIFF_FIRST_ORDER_FUNCTION_H_
+#define CERES_PUBLIC_NUMERIC_DIFF_FIRST_ORDER_FUNCTION_H_
+
+#include <algorithm>
+#include <memory>
+
+#include "ceres/first_order_function.h"
+#include "ceres/internal/eigen.h"
+#include "ceres/internal/fixed_array.h"
+#include "ceres/internal/numeric_diff.h"
+#include "ceres/internal/parameter_dims.h"
+#include "ceres/internal/variadic_evaluate.h"
+#include "ceres/numeric_diff_options.h"
+#include "ceres/types.h"
+
+
+namespace ceres {
+
+// Creates FirstOrderFunctions as needed by the GradientProblem
+// framework, with gradients computed via numeric differentiation. For
+// more information on numeric differentiation, see the wikipedia
+// article at https://en.wikipedia.org/wiki/Numerical_differentiation
+//
+// To get an numerically differentiated cost function, you must define
+// a class with an operator() (a functor) that computes the cost.
+//
+// The function must write the computed value in the last argument
+// (the only non-const one) and return true to indicate success.
+//
+// For example, consider a scalar error e = x'y - a, where both x and y are
+// two-dimensional column vector parameters, the prime sign indicates
+// transposition, and a is a constant.
+//
+// To write an numerically-differentiable cost function for the above model,
+// first define the object
+//
+//  class QuadraticCostFunctor {
+//   public:
+//    explicit QuadraticCostFunctor(double a) : a_(a) {}
+//    bool operator()(const double* const xy, double* cost) const {
+//      constexpr int kInputVectorLength = 2;
+//      const double* const x = xy;
+//      const double* const y = xy + kInputVectorLength;
+//      *cost = x[0] * y[0] + x[1] * y[1] - a_;
+//      return true;
+//    }
+//
+//   private:
+//    double a_;
+//  };
+//
+//
+// Note that in the declaration of operator() the input parameters xy
+// come first, and are passed as const pointers to array of
+// doubles. The output cost is the last parameter.
+//
+// Then given this class definition, the numerically differentiated
+// first order function with central differences used for computing the
+// derivative can be constructed as follows.
+//
+//   FirstOrderFunction* function
+//       = new NumericDiffFirstOrderFunction<MyScalarCostFunctor, CENTRAL, 4>(
+//           new QuadraticCostFunctor(1.0));                   ^     ^     ^
+//                                                             |     |     |
+//                                 Finite Differencing Scheme -+     |     |
+//                                 Dimension of xy ------------------------+
+//
+//
+// In the instantiation above, the template parameters following
+// "QuadraticCostFunctor", "CENTRAL, 4", describe the finite
+// differencing scheme as "central differencing" and the functor as
+// computing its cost from a 4 dimensional input.
+template <typename FirstOrderFunctor,
+          NumericDiffMethodType method,
+          int kNumParameters>
+class NumericDiffFirstOrderFunction : public FirstOrderFunction {
+ public:
+  NumericDiffFirstOrderFunction(
+      FirstOrderFunctor* functor,
+      Ownership ownership = TAKE_OWNERSHIP,
+      const NumericDiffOptions& options = NumericDiffOptions())
+      : functor_(functor), ownership_(ownership), options_(options) {
+    static_assert(kNumParameters > 0, "kNumParameters must be positive");
+  }
+
+  ~NumericDiffFirstOrderFunction() override {
+    if (ownership_ != TAKE_OWNERSHIP) {
+      functor_.release();
+    }
+  }
+
+  bool Evaluate(const double* const parameters,
+                double* cost,
+                double* gradient) const override {
+    using ParameterDims = internal::StaticParameterDims<kNumParameters>;
+    constexpr int kNumResiduals = 1;
+
+    // Get the function value (cost) at the the point to evaluate.
+    if (!internal::VariadicEvaluate<ParameterDims>(
+            *functor_, &parameters, cost)) {
+      return false;
+    }
+
+    if (gradient == nullptr) {
+      return true;
+    }
+
+    // Create a copy of the parameters which will get mutated.
+    internal::FixedArray<double, 32> parameters_copy(kNumParameters);
+    std::copy_n(parameters, kNumParameters, parameters_copy.data());
+    double* parameters_ptr = parameters_copy.data();
+    internal::EvaluateJacobianForParameterBlocks<
+        ParameterDims>::template Apply<method, kNumResiduals>(functor_.get(),
+                                                              cost,
+                                                              options_,
+                                                              kNumResiduals,
+                                                              &parameters_ptr,
+                                                              &gradient);
+    return true;
+  }
+
+  int NumParameters() const override { return kNumParameters; }
+
+ private:
+  std::unique_ptr<FirstOrderFunctor> functor_;
+  Ownership ownership_;
+  NumericDiffOptions options_;
+};
+
+}  // namespace ceres
+
+#endif  // CERES_PUBLIC_NUMERIC_DIFF_FIRST_ORDER_FUNCTION_H_
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 80f8bdc..18f6090 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -486,6 +486,7 @@
   ceres_test(minimizer)
   ceres_test(normal_prior)
   ceres_test(numeric_diff_cost_function)
+  ceres_test(numeric_diff_first_order_function)
   ceres_test(ordered_groups)
   ceres_test(parallel_for)
   ceres_test(parallel_utils)
@@ -549,4 +550,3 @@
 
   add_subdirectory(autodiff_benchmarks)
 endif (BUILD_BENCHMARKS)
-
diff --git a/internal/ceres/autodiff_first_order_function_test.cc b/internal/ceres/autodiff_first_order_function_test.cc
index 7db7835..4870a45 100644
--- a/internal/ceres/autodiff_first_order_function_test.cc
+++ b/internal/ceres/autodiff_first_order_function_test.cc
@@ -44,7 +44,7 @@
   explicit QuadraticCostFunctor(double a) : a_(a) {}
   template <typename T>
   bool operator()(const T* const x, T* cost) const {
-    cost[0] = x[0] * x[1] + x[2] * x[3] - T(a_);
+    cost[0] = x[0] * x[1] + x[2] * x[3] - a_;
     return true;
   }
 
diff --git a/internal/ceres/numeric_diff_first_order_function_test.cc b/internal/ceres/numeric_diff_first_order_function_test.cc
new file mode 100644
index 0000000..1926caf
--- /dev/null
+++ b/internal/ceres/numeric_diff_first_order_function_test.cc
@@ -0,0 +1,79 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2021 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/numeric_diff_first_order_function.h"
+
+#include <memory>
+
+#include "ceres/array_utils.h"
+#include "ceres/first_order_function.h"
+#include "gtest/gtest.h"
+
+namespace ceres {
+namespace internal {
+
+class QuadraticCostFunctor {
+ public:
+  explicit QuadraticCostFunctor(double a) : a_(a) {}
+  bool operator()(const double* const x, double* cost) const {
+    cost[0] = x[0] * x[1] + x[2] * x[3] - a_;
+    return true;
+  }
+
+ private:
+  double a_;
+};
+
+TEST(NumericDiffFirstOrderFunction, BilinearDifferentiationTest) {
+  std::unique_ptr<FirstOrderFunction> function(
+      new NumericDiffFirstOrderFunction<QuadraticCostFunctor, CENTRAL, 4>(
+          new QuadraticCostFunctor(1.0)));
+
+  double parameters[4] = {1.0, 2.0, 3.0, 4.0};
+  double gradient[4];
+  double cost;
+
+  function->Evaluate(parameters, &cost, nullptr);
+  EXPECT_EQ(cost, 13.0);
+
+  cost = -1.0;
+  function->Evaluate(parameters, &cost, gradient);
+
+  EXPECT_EQ(cost, 13.0);
+
+  const double kTolerance = 1e-9;
+  EXPECT_NEAR(gradient[0], parameters[1], kTolerance);
+  EXPECT_NEAR(gradient[1], parameters[0], kTolerance);
+  EXPECT_NEAR(gradient[2], parameters[3], kTolerance);
+  EXPECT_NEAR(gradient[3], parameters[2], kTolerance);
+}
+
+}  // namespace internal
+}  // namespace ceres