diff --git a/math/mathcore/inc/Fit/Chi2FCN.h b/math/mathcore/inc/Fit/Chi2FCN.h
index bc72f64d680c5..f45cb872052ea 100644
--- a/math/mathcore/inc/Fit/Chi2FCN.h
+++ b/math/mathcore/inc/Fit/Chi2FCN.h
@@ -126,6 +126,14 @@ class Chi2FCN : public BasicFCN<DerivFunType, ModelFunType, BinData> {
       FitUtil::Evaluate<T>::EvalChi2Gradient(BaseFCN::ModelFunction(), BaseFCN::Data(), x, g, fNEffPoints,
                                              fExecutionPolicy);
    }
+   /// In some cases, the gradient algorithm will use information from the previous step, these can be passed
+   /// in with this overload. The `previous_*` arrays can also be used to return second derivative and step size
+   /// so that these can be passed forward again as well at the call site, if necessary.
+   /// \warning This implementation just calls the two-parameter overload.
+   virtual void Gradient(const double *x, double *g, double */*previous_grad*/, double */*previous_g2*/, double */*previous_gstep*/) const
+   {
+      Gradient(x, g);
+   }
 
    /// get type of fit method function
    virtual  typename BaseObjFunction::Type_t Type() const { return BaseObjFunction::kLeastSquare; }
@@ -154,6 +162,15 @@ class Chi2FCN : public BasicFCN<DerivFunType, ModelFunType, BinData> {
       Gradient(x, fGrad.data());
       return fGrad[icoord];
    }
+   /// In some cases, the derivative algorithm will use information from the previous step, these can be passed
+   /// in with this overload. The `previous_*` arrays can also be used to return second derivative and step size
+   /// so that these can be passed forward again as well at the call site, if necessary.
+   /// \warning This implementation just calls the two-parameter overload.
+   virtual double DoDerivative(const double *x, unsigned int icoord, double * /*previous_grad*/, double * /*previous_g2*/,
+                               double * /*previous_gstep*/) const
+   {
+      return DoDerivative(x, icoord);
+   }
 
 
    mutable unsigned int fNEffPoints;  // number of effective points used in the fit
diff --git a/math/mathcore/inc/Fit/LogLikelihoodFCN.h b/math/mathcore/inc/Fit/LogLikelihoodFCN.h
index 02b9ec1715a2b..bada8b1d49225 100644
--- a/math/mathcore/inc/Fit/LogLikelihoodFCN.h
+++ b/math/mathcore/inc/Fit/LogLikelihoodFCN.h
@@ -129,6 +129,14 @@ class LogLikelihoodFCN : public BasicFCN<DerivFunType,ModelFunType,UnBinData>  {
       FitUtil::Evaluate<typename BaseFCN::T>::EvalLogLGradient(BaseFCN::ModelFunction(), BaseFCN::Data(), x, g,
                                                                fNEffPoints, fExecutionPolicy);
    }
+   /// In some cases, the gradient algorithm will use information from the previous step, these can be passed
+   /// in with this overload. The `previous_*` arrays can also be used to return second derivative and step size
+   /// so that these can be passed forward again as well at the call site, if necessary.
+   /// \warning This overload just calls the two-parameter version.
+   virtual void Gradient(const double *x, double *g, double */*previous_grad*/, double */*previous_g2*/, double */*previous_gstep*/) const
+   {
+      Gradient(x, g);
+   }
 
    /// get type of fit method function
    virtual  typename BaseObjFunction::Type_t Type() const { return BaseObjFunction::kLogLikelihood; }
@@ -162,6 +170,14 @@ class LogLikelihoodFCN : public BasicFCN<DerivFunType,ModelFunType,UnBinData>  {
       Gradient(x, &fGrad[0]);
       return fGrad[icoord];
    }
+   /// In some cases, the derivative algorithm will use information from the previous step, these can be passed
+   /// in with this overload. The `previous_*` arrays can also be used to return second derivative and step size
+   /// so that these can be passed forward again as well at the call site, if necessary.
+   virtual double DoDerivative(const double *x, unsigned int icoord, double * /*previous_grad*/, double * /*previous_g2*/,
+                               double * /*previous_gstep*/) const
+   {
+      return DoDerivative(x, icoord);
+   }
 
 
       //data member
diff --git a/math/mathcore/inc/Fit/PoissonLikelihoodFCN.h b/math/mathcore/inc/Fit/PoissonLikelihoodFCN.h
index 85d1025b0f6a9..7cc2ca3aeb71c 100644
--- a/math/mathcore/inc/Fit/PoissonLikelihoodFCN.h
+++ b/math/mathcore/inc/Fit/PoissonLikelihoodFCN.h
@@ -130,6 +130,14 @@ class PoissonLikelihoodFCN : public BasicFCN<DerivFunType,ModelFunType,BinData>
       FitUtil::Evaluate<typename BaseFCN::T>::EvalPoissonLogLGradient(BaseFCN::ModelFunction(), BaseFCN::Data(), x, g,
                                                                       fNEffPoints, fExecutionPolicy);
    }
+   /// In some cases, the gradient algorithm will use information from the previous step, these can be passed
+   /// in with this overload. The `previous_*` arrays can also be used to return second derivative and step size
+   /// so that these can be passed forward again as well at the call site, if necessary.
+   /// \warning This implementation just calls the two-parameter overload.
+   virtual void Gradient(const double *x, double *g, double */*previous_grad*/, double */*previous_g2*/, double */*previous_gstep*/) const
+   {
+      Gradient(x, g);
+   }
 
    /// get type of fit method function
    virtual  typename BaseObjFunction::Type_t Type() const { return BaseObjFunction::kLogLikelihood; }
@@ -170,6 +178,15 @@ class PoissonLikelihoodFCN : public BasicFCN<DerivFunType,ModelFunType,BinData>
       Gradient(x, &fGrad[0]);
       return fGrad[icoord];
    }
+   /// In some cases, the derivative algorithm will use information from the previous step, these can be passed
+   /// in with this overload. The `previous_*` arrays can also be used to return second derivative and step size
+   /// so that these can be passed forward again as well at the call site, if necessary.
+   /// \warning This implementation just calls the two-parameter overload.
+   virtual double DoDerivative(const double *x, unsigned int icoord, double * /*previous_grad*/, double * /*previous_g2*/,
+                               double * /*previous_gstep*/) const
+   {
+      return DoDerivative(x, icoord);
+   }
 
 
       //data member
diff --git a/math/mathcore/inc/Math/Functor.h b/math/mathcore/inc/Math/Functor.h
index 43214aa8d2190..0ac2fbc9605a6 100644
--- a/math/mathcore/inc/Math/Functor.h
+++ b/math/mathcore/inc/Math/Functor.h
@@ -128,6 +128,14 @@ private :
       return fFunc.Derivative(x,icoord);
    }
 
+   // TODO: implementing this will require extending with another function signature in GradFunctor
+   /// \warning This overload just calls the two-parameter version.
+   inline double DoDerivative(const double *x, unsigned int icoord, double */*previous_grad*/, double */*previous_g2*/,
+                              double */*previous_gstep*/) const
+   {
+      return DoDerivative(x, icoord);
+   }
+
 
    unsigned int fDim;
    mutable Func fFunc;  // should here be a reference and pass a non-const ref in ctor
@@ -211,6 +219,13 @@ private :
       return fGradFunc(x, icoord);
    }
 
+   // TODO: implementing this will require extending with another function signature in GradFunctor
+   /// \warning This overload just calls the two-parameter version.
+   inline double DoDerivative(const double *x, unsigned int icoord, double */*previous_grad*/, double */*previous_g2*/,
+                              double */*previous_gstep*/) const
+   {
+      return DoDerivative(x, icoord);
+   }
 
    unsigned int fDim;
    mutable Func fFunc;
@@ -348,6 +363,14 @@ private :
       return ((*fObj).*fGradMemFn)(x,icoord);
    }
 
+   // TODO: implementing this will require extending with another function signature in GradFunctor
+   /// \warning This overload just calls the two-parameter version.
+   inline double DoDerivative(const double *x, unsigned int icoord, double */*previous_grad*/, double */*previous_g2*/,
+                              double */*previous_gstep*/) const
+   {
+      return DoDerivative(x, icoord);
+   }
+
    unsigned int fDim;
    mutable PointerToObj fObj;
    PointerToMemFn fMemFn;
@@ -702,6 +725,14 @@ private :
       return fImpl->Derivative(x,icoord);
    }
 
+   // TODO: implementing this will require extending with another function signature in GradFunctor
+   /// \warning This overload just calls the two-parameter version.
+   inline double DoDerivative(const double *x, unsigned int icoord, double */*previous_grad*/, double */*previous_g2*/,
+                              double */*previous_gstep*/) const
+   {
+      return DoDerivative(x, icoord);
+   }
+
    std::unique_ptr<Impl> fImpl;    // pointer to base grad functor handler
 
 
diff --git a/math/mathcore/inc/Math/IFunction.h b/math/mathcore/inc/Math/IFunction.h
index e7871c4f06f76..474a5aed6137d 100644
--- a/math/mathcore/inc/Math/IFunction.h
+++ b/math/mathcore/inc/Math/IFunction.h
@@ -215,6 +215,14 @@ namespace ROOT {
             Return the partial derivative with respect to the passed coordinate
          */
          T Derivative(const T *x, unsigned int icoord = 0) const { return DoDerivative(x, icoord); }
+         /// In some cases, the derivative algorithm will use information from the previous step, these can be passed
+         /// in with this overload. The `previous_*` arrays can also be used to return second derivative and step size
+         /// so that these can be passed forward again as well at the call site, if necessary.
+         T Derivative(const T *x, unsigned int icoord, T *previous_grad, T *previous_g2,
+                      T *previous_gstep) const
+         {
+            return DoDerivative(x, icoord, previous_grad, previous_g2, previous_gstep);
+         }
 
          /**
              Optimized method to evaluate at the same time the function value and derivative at a point x.
@@ -232,6 +240,14 @@ namespace ROOT {
             function to evaluate the derivative with respect each coordinate. To be implemented by the derived class
          */
          virtual T DoDerivative(const T *x, unsigned int icoord) const = 0;
+         /// In some cases, the derivative algorithm will use information from the previous step, these can be passed
+         /// in with this overload. The `previous_*` arrays can also be used to return second derivative and step size
+         /// so that these can be passed forward again as well at the call site, if necessary.
+         virtual T DoDerivative(const T *x, unsigned int icoord, T * /*previous_grad*/, T * /*previous_g2*/,
+                                T * /*previous_gstep*/) const
+         {
+            return DoDerivative(x, icoord);
+         };
       };
 
 //___________________________________________________________________________________
@@ -346,6 +362,16 @@ namespace ROOT {
                grad[icoord] = BaseGrad::Derivative(x, icoord);
          }
 
+         /// In some cases, the gradient algorithm will use information from the previous step, these can be passed
+         /// in with this overload. The `previous_*` arrays can also be used to return second derivative and step size
+         /// so that these can be passed forward again as well at the call site, if necessary.
+         virtual void Gradient(const T *x, T *grad, T *previous_grad, T *previous_g2, T *previous_gstep) const
+         {
+            unsigned int ndim = NDim();
+            for (unsigned int icoord  = 0; icoord < ndim; ++icoord)
+               grad[icoord] = BaseGrad::Derivative(x, icoord, previous_grad, previous_g2, previous_gstep);
+         }
+
          using  BaseFunc::NDim;
 
          /**
diff --git a/math/mathcore/inc/Math/MinimTransformFunction.h b/math/mathcore/inc/Math/MinimTransformFunction.h
index 94d1877778326..c8eac81b3e4d7 100644
--- a/math/mathcore/inc/Math/MinimTransformFunction.h
+++ b/math/mathcore/inc/Math/MinimTransformFunction.h
@@ -123,6 +123,14 @@ class MinimTransformFunction : public IMultiGradFunction {
       //std::cout << "Derivative icoord (ext)" << fIndex[icoord] << "   dtrafo " << dExtdInt << "  " << deriv << std::endl;
       return deriv * dExtdInt;
    }
+   /// In some cases, the derivative algorithm will use information from the previous step, these can be passed
+   /// in with this overload. The `previous_*` arrays can also be used to return second derivative and step size
+   /// so that these can be passed forward again as well at the call site, if necessary.
+   virtual double DoDerivative(const double *x, unsigned int icoord, double * /*previous_grad*/, double * /*previous_g2*/,
+                               double * /*previous_gstep*/) const
+   {
+      return DoDerivative(x, icoord);
+   }
 
    // copy constructor for this class (disable by having it private)
    MinimTransformFunction( const MinimTransformFunction & ) :
diff --git a/math/mathcore/test/fit/testMinim.cxx b/math/mathcore/test/fit/testMinim.cxx
index 27afa917ab1a1..f0275e28eb486 100644
--- a/math/mathcore/test/fit/testMinim.cxx
+++ b/math/mathcore/test/fit/testMinim.cxx
@@ -188,6 +188,10 @@ public :
       }
 
    }
+   void Gradient(const double *x, double *g, double */*previous_grad*/, double */*previous_g2*/, double */*previous_gstep*/) const
+   {
+      Gradient(x, g);
+   }
 
 #ifdef USE_FDF
    void FdF (const double * x, double & f, double * g) const {
@@ -240,6 +244,11 @@ public :
       Gradient(x,&g[0]);
       return  g[i];
    }
+   double DoDerivative(const double *x, unsigned int icoord, double * /*previous_grad*/, double * /*previous_g2*/,
+                       double * /*previous_gstep*/) const
+   {
+      return DoDerivative(x, icoord);
+   }
 
 private:
 
@@ -322,6 +331,10 @@ public :
 
 
    }
+   void Gradient(const double *x, double *g, double */*previous_grad*/, double */*previous_g2*/, double */*previous_gstep*/) const
+   {
+      Gradient(x, g);
+   }
 
    private:
 
@@ -342,6 +355,11 @@ public :
       Gradient(x,&g[0]);
       return  g[i];
    }
+   double DoDerivative(const double *x, unsigned int icoord, double * /*previous_grad*/, double * /*previous_g2*/,
+                       double * /*previous_gstep*/) const
+   {
+      return DoDerivative(x, icoord);
+   }
 
    void DoCalculatefi(const double * x) const {
       // calculate the i- element ; F(X) = Sum {fi]
diff --git a/math/mathmore/inc/Math/GSLNLSMinimizer.h b/math/mathmore/inc/Math/GSLNLSMinimizer.h
index 46cd3364c5bb2..8bc2a9d0f311a 100644
--- a/math/mathmore/inc/Math/GSLNLSMinimizer.h
+++ b/math/mathmore/inc/Math/GSLNLSMinimizer.h
@@ -102,6 +102,14 @@ class LSResidualFunc : public IMultiGradFunction {
       double f0 = 0;
       FdF(x,f0,g);
    }
+   /// In some cases, the gradient algorithm will use information from the previous step, these can be passed
+   /// in with this overload. The `previous_*` arrays can also be used to return second derivative and step size
+   /// so that these can be passed forward again as well at the call site, if necessary.
+   /// \warning This implementation just calls the two-parameter overload.
+   virtual void Gradient(const double *x, double *g, double */*previous_grad*/, double */*previous_g2*/, double */*previous_gstep*/) const
+   {
+      Gradient(x, g);
+   }
 
    void FdF (const double * x, double & f, double * g) const {
       unsigned int n = NDim();
@@ -129,6 +137,15 @@ class LSResidualFunc : public IMultiGradFunction {
       fX2[icoord] += kEps;
       return ( DoEval(&fX2.front()) - DoEval(x) )/kEps;
    }
+   /// In some cases, the derivative algorithm will use information from the previous step, these can be passed
+   /// in with this overload. The `previous_*` arrays can also be used to return second derivative and step size
+   /// so that these can be passed forward again as well at the call site, if necessary.
+   /// \warning This implementation just calls the two-parameter overload.
+   virtual double DoDerivative(const double *x, unsigned int icoord, double * /*previous_grad*/, double * /*previous_g2*/,
+                               double * /*previous_gstep*/) const
+   {
+      return DoDerivative(x, icoord);
+   }
 
    unsigned int fIndex;
    const ROOT::Math::FitMethodFunction * fChi2;
diff --git a/math/mathmore/inc/Math/MultiNumGradFunction.h b/math/mathmore/inc/Math/MultiNumGradFunction.h
index 2bec540179222..1676661fdbde8 100644
--- a/math/mathmore/inc/Math/MultiNumGradFunction.h
+++ b/math/mathmore/inc/Math/MultiNumGradFunction.h
@@ -122,6 +122,15 @@ class MultiNumGradFunction : public IMultiGradFunction {
 
    // calculate derivative using mathcore derivator
    double DoDerivative (const double * x, unsigned int icoord  ) const;
+   /// In some cases, the derivative algorithm will use information from the previous step, these can be passed
+   /// in with this overload. The `previous_*` arrays can also be used to return second derivative and step size
+   /// so that these can be passed forward again as well at the call site, if necessary.
+   /// \warning This implementation just calls the two-parameter overload.
+   virtual double DoDerivative(const double *x, unsigned int icoord, double * /*previous_grad*/, double * /*previous_g2*/,
+                               double * /*previous_gstep*/) const
+   {
+      return DoDerivative(x, icoord);
+   }
 
    // adapat internal function type to IMultiGenFunction needed by derivative calculation
    const IMultiGenFunction * fFunc;
diff --git a/math/minuit2/inc/Minuit2/FCNGradAdapter.h b/math/minuit2/inc/Minuit2/FCNGradAdapter.h
index 6d42bab48f340..5e186918da0c9 100644
--- a/math/minuit2/inc/Minuit2/FCNGradAdapter.h
+++ b/math/minuit2/inc/Minuit2/FCNGradAdapter.h
@@ -54,6 +54,18 @@ class FCNGradAdapter : public FCNGradientBase {
       });
       return fGrad;
    }
+   std::vector<double> Gradient(const std::vector<double> &v, double *previous_grad, double *previous_g2,
+                                double *previous_gstep) const override
+   {
+      fFunc.Gradient(&v[0], &fGrad[0], previous_grad, previous_g2, previous_gstep);
+
+      MnPrint("FCNGradAdapter").Debug([&](std::ostream &os) {
+         os << "gradient in FCNAdapter = {";
+         for (unsigned int i = 0; i < fGrad.size(); ++i)
+            os << fGrad[i] << (i == fGrad.size() - 1 ? '}' : '\t');
+      });
+      return fGrad;
+   }
    // forward interface
    // virtual double operator()(int npar, double* params,int iflag = 4) const;
    bool CheckGradient() const override { return false; }
diff --git a/math/minuit2/inc/Minuit2/FCNGradientBase.h b/math/minuit2/inc/Minuit2/FCNGradientBase.h
index 4300528b65929..cc7259fd57dbc 100644
--- a/math/minuit2/inc/Minuit2/FCNGradientBase.h
+++ b/math/minuit2/inc/Minuit2/FCNGradientBase.h
@@ -41,6 +41,8 @@ class FCNGradientBase : public FCNBase {
    virtual ~FCNGradientBase() {}
 
    virtual std::vector<double> Gradient(const std::vector<double> &) const = 0;
+   virtual std::vector<double> Gradient(const std::vector<double> &parameters, double */*previous_grad*/, double */*previous_g2*/,
+                                        double */*previous_gstep*/) const { return Gradient(parameters); };
 
    virtual bool CheckGradient() const { return true; }
 
diff --git a/math/minuit2/inc/Minuit2/NumericalDerivator.h b/math/minuit2/inc/Minuit2/NumericalDerivator.h
index ee6203274f066..e2c90731b68c2 100644
--- a/math/minuit2/inc/Minuit2/NumericalDerivator.h
+++ b/math/minuit2/inc/Minuit2/NumericalDerivator.h
@@ -70,8 +70,7 @@ class NumericalDerivator {
    double Ext2int(const ROOT::Fit::ParameterSettings &parameter, double val) const;
    double DInt2Ext(const ROOT::Fit::ParameterSettings &parameter, double val) const;
 
-   void SetInitialGradient(const ROOT::Math::IBaseFunctionMultiDim *function,
-                           const std::vector<ROOT::Fit::ParameterSettings> &parameters,
+   void SetInitialGradient(const std::vector<ROOT::Fit::ParameterSettings> &parameters,
                            std::vector<DerivatorElement> &gradient);
 
    inline bool AlwaysExactlyMimicMinuit2() const { return fAlwaysExactlyMimicMinuit2; }
diff --git a/math/minuit2/src/ExternalInternalGradientCalculator.cxx b/math/minuit2/src/ExternalInternalGradientCalculator.cxx
index 19700394fc9d7..0f93f2bc6e5e8 100644
--- a/math/minuit2/src/ExternalInternalGradientCalculator.cxx
+++ b/math/minuit2/src/ExternalInternalGradientCalculator.cxx
@@ -13,13 +13,13 @@
 #include "Minuit2/MnUserTransformation.h"
 #include "Minuit2/FunctionGradient.h"
 #include "Minuit2/MinimumParameters.h"
+#include "Minuit2/MnPrint.h"
 
 namespace ROOT {
 namespace Minuit2 {
 
 FunctionGradient ExternalInternalGradientCalculator::operator()(const MinimumParameters &par) const
 {
-   // evaluate analytical gradient. take care of parameter transformations
    std::vector<double> par_vec;
    par_vec.resize(par.Vec().size());
    for (std::size_t ix = 0; ix < par.Vec().size(); ++ix) {
@@ -35,18 +35,42 @@ FunctionGradient ExternalInternalGradientCalculator::operator()(const MinimumPar
       v(i) = grad[ext];
    }
 
-#ifdef DEBUG
-   std::cout << "User given gradient in Minuit2" << v << std::endl;
-#endif
+   MnPrint print("ExternalInternalGradientCalculator");
+   print.Debug("User given gradient in Minuit2", v);
 
    return FunctionGradient(v);
 }
 
 FunctionGradient
-ExternalInternalGradientCalculator::operator()(const MinimumParameters &par, const FunctionGradient &) const
+ExternalInternalGradientCalculator::operator()(const MinimumParameters &par, const FunctionGradient &functionGradient) const
 {
-   // needed from base class
-   return (*this)(par);
+   std::vector<double> par_vec;
+   par_vec.resize(par.Vec().size());
+   for (std::size_t ix = 0; ix < par.Vec().size(); ++ix) {
+      par_vec[ix] = par.Vec()(ix);
+   }
+
+   std::vector<double> previous_grad(functionGradient.Grad().Data(), functionGradient.Grad().Data() + functionGradient.Grad().size());
+   std::vector<double> previous_g2(functionGradient.G2().Data(), functionGradient.G2().Data() + functionGradient.G2().size());
+   std::vector<double> previous_gstep(functionGradient.Gstep().Data(), functionGradient.Gstep().Data() + functionGradient.Gstep().size());
+
+   std::vector<double> grad = fGradCalc.Gradient(par_vec, previous_grad.data(), previous_g2.data(), previous_gstep.data());
+   assert(grad.size() == fTransformation.Parameters().size());
+
+   MnAlgebraicVector v(par.Vec().size());
+   MnAlgebraicVector v_g2(par.Vec().size());
+   MnAlgebraicVector v_gstep(par.Vec().size());
+   for (unsigned int i = 0; i < par.Vec().size(); i++) {
+      unsigned int ext = fTransformation.ExtOfInt(i);
+      v(i) = grad[ext];
+      v_g2(i) = previous_g2[ext];
+      v_gstep(i) = previous_gstep[ext];
+   }
+
+   MnPrint print("ExternalInternalGradientCalculator");
+   print.Debug("User given gradient in Minuit2", v, "g2", v_g2, "step size", v_gstep);
+
+   return FunctionGradient(v, v_g2, v_gstep);
 }
 
 } // namespace Minuit2
diff --git a/math/minuit2/src/NumericalDerivator.cxx b/math/minuit2/src/NumericalDerivator.cxx
index 618334d0e55a4..758cd6c8d8c28 100644
--- a/math/minuit2/src/NumericalDerivator.cxx
+++ b/math/minuit2/src/NumericalDerivator.cxx
@@ -221,8 +221,7 @@ double NumericalDerivator::DInt2Ext(const ROOT::Fit::ParameterSettings &paramete
 // MODIFIED:
 /// This function was not implemented as in Minuit2. Now it copies the behavior
 /// of InitialGradientCalculator. See https://github.com/roofit-dev/root/issues/10
-void NumericalDerivator::SetInitialGradient(const ROOT::Math::IBaseFunctionMultiDim *,
-                                            const std::vector<ROOT::Fit::ParameterSettings> &parameters,
+void NumericalDerivator::SetInitialGradient(const std::vector<ROOT::Fit::ParameterSettings> &parameters,
                                             std::vector<DerivatorElement> &gradient)
 {
    // set an initial gradient using some given steps
diff --git a/math/minuit2/test/MnTutorial/Quad1F.h b/math/minuit2/test/MnTutorial/Quad1F.h
index ba6fabf1a93a0..276282e7ec004 100644
--- a/math/minuit2/test/MnTutorial/Quad1F.h
+++ b/math/minuit2/test/MnTutorial/Quad1F.h
@@ -37,6 +37,8 @@ class Quad1F : public FCNGradientBase {
 
       return (std::vector<double>(1, 2. * x));
    }
+   virtual std::vector<double> Gradient(const std::vector<double> &parameters, double */*previous_grad*/, double */*previous_g2*/,
+                                        double */*previous_gstep*/) const { return Gradient(parameters); };
 
    void SetErrorDef(double up) { fErrorDef = up; }
 
diff --git a/math/minuit2/test/MnTutorial/Quad4F.h b/math/minuit2/test/MnTutorial/Quad4F.h
index 1089184ca7a69..ed90fc4925079 100644
--- a/math/minuit2/test/MnTutorial/Quad4F.h
+++ b/math/minuit2/test/MnTutorial/Quad4F.h
@@ -72,6 +72,8 @@ class Quad4FGrad : public FCNGradientBase {
       g[3] = 2. * w;
       return g;
    }
+   virtual std::vector<double> Gradient(const std::vector<double> &parameters, double */*previous_grad*/, double */*previous_g2*/,
+                                        double */*previous_gstep*/) const { return Gradient(parameters); };
 
    double Up() const { return 1.; }
 
diff --git a/math/minuit2/test/testMinimizer.cxx b/math/minuit2/test/testMinimizer.cxx
index 92f7817204b1a..dea116e13fd85 100644
--- a/math/minuit2/test/testMinimizer.cxx
+++ b/math/minuit2/test/testMinimizer.cxx
@@ -169,6 +169,10 @@ class TrigoFletcherFunction : public ROOT::Math::IMultiGradFunction {
          }
       }
    }
+   void Gradient(const double *x, double *g, double */*previous_grad*/, double */*previous_g2*/, double */*previous_gstep*/) const
+   {
+      Gradient(x, g);
+   }
 
 #ifdef USE_FDF
    void FdF(const double *x, double &f, double *g) const
@@ -220,6 +224,11 @@ class TrigoFletcherFunction : public ROOT::Math::IMultiGradFunction {
       Gradient(x, &g[0]);
       return g[i];
    }
+   double DoDerivative(const double *x, unsigned int icoord, double * /*previous_grad*/, double * /*previous_g2*/,
+                       double * /*previous_gstep*/) const
+   {
+      return DoDerivative(x, icoord);
+   }
 
 private:
    unsigned int fDim;
@@ -286,6 +295,10 @@ class ChebyQuadFunction : public ROOT::Math::IMultiGradFunction {
          g[j] = 2. * g[j] / double(n);
       }
    }
+   void Gradient(const double *x, double *g, double */*previous_grad*/, double */*previous_g2*/, double */*previous_gstep*/) const
+   {
+      Gradient(x, g);
+   }
 
 private:
    double DoEval(const double *x) const
@@ -306,6 +319,11 @@ class ChebyQuadFunction : public ROOT::Math::IMultiGradFunction {
       Gradient(x, &g[0]);
       return g[i];
    }
+   double DoDerivative(const double *x, unsigned int icoord, double * /*previous_grad*/, double * /*previous_g2*/,
+                       double * /*previous_gstep*/) const
+   {
+      return DoDerivative(x, icoord);
+   }
 
    void DoCalculatefi(const double *x) const
    {
diff --git a/roofit/roofitcore/inc/RooGradMinimizerFcn.h b/roofit/roofitcore/inc/RooGradMinimizerFcn.h
index 829fc228d45d5..f17563b60e0f6 100644
--- a/roofit/roofitcore/inc/RooGradMinimizerFcn.h
+++ b/roofit/roofitcore/inc/RooGradMinimizerFcn.h
@@ -70,6 +70,8 @@ class RooGradMinimizerFcn : public ROOT::Math::IMultiGradFunction, public RooAbs
    // IMultiGradFunction overrides
    double DoEval(const double *x) const override;
    double DoDerivative(const double *x, unsigned int icoord) const override;
+   double DoDerivative(const double *x, unsigned int i_component, double *previous_grad,
+                       double *previous_g2, double *previous_gstep) const override;
 
    // members
    // mutable because ROOT::Math::IMultiGradFunction::DoDerivative is const
diff --git a/roofit/roofitcore/src/RooGradMinimizerFcn.cxx b/roofit/roofitcore/src/RooGradMinimizerFcn.cxx
index 55e5d66ab5010..6e389a1d555ac 100644
--- a/roofit/roofitcore/src/RooGradMinimizerFcn.cxx
+++ b/roofit/roofitcore/src/RooGradMinimizerFcn.cxx
@@ -64,7 +64,7 @@ ROOT::Math::IMultiGradFunction *RooGradMinimizerFcn::Clone() const
 void RooGradMinimizerFcn::synchronizeGradientParameterSettings(
    std::vector<ROOT::Fit::ParameterSettings> &parameter_settings) const
 {
-   _gradf.SetInitialGradient(_context->getMultiGenFcn(), parameter_settings, _grad);
+   _gradf.SetInitialGradient(parameter_settings, _grad);
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -215,6 +215,18 @@ double RooGradMinimizerFcn::DoDerivative(const double *x, unsigned int i_compone
    return _grad[i_component].derivative;
 }
 
+double RooGradMinimizerFcn::DoDerivative(const double *x, unsigned int i_component, double *previous_grad,
+                                         double *previous_g2, double *previous_gstep) const
+{
+   syncParameters(x);
+   _grad[i_component] = {previous_grad[i_component], previous_g2[i_component], previous_gstep[i_component]};
+   runDerivator(i_component);
+   previous_grad[i_component] = _grad[i_component].derivative;
+   previous_g2[i_component] = _grad[i_component].second_derivative;
+   previous_gstep[i_component] = _grad[i_component].step_size;
+   return _grad[i_component].derivative;
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 
 void RooGradMinimizerFcn::setStrategy(int istrat)
diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt
index 2ebcb7507e091..890723e81aefa 100644
--- a/roofit/roofitcore/test/CMakeLists.txt
+++ b/roofit/roofitcore/test/CMakeLists.txt
@@ -42,6 +42,5 @@ if(NOT MSVC OR win_broken_tests)
 endif()
 ROOT_ADD_GTEST(testRooProductPdf testRooProductPdf.cxx LIBRARIES RooFitCore)
 ROOT_ADD_GTEST(testNaNPacker testNaNPacker.cxx LIBRARIES RooFitCore)
-ROOT_ADD_GTEST(testRooSimultaneous testRooSimultaneous.cxx LIBRARIES RooFitCore)
-#ROOT_ADD_GTEST(testRooGradMinimizerFcn testRooGradMinimizerFcn.cxx LIBRARIES RooFitCore)
-
+ROOT_ADD_GTEST(testRooSimultaneous testRooSimultaneous.cxx LIBRARIES RooFitCore RooFit)
+ROOT_ADD_GTEST(testRooGradMinimizerFcn testRooGradMinimizerFcn.cxx LIBRARIES RooFitCore)
diff --git a/roofit/roofitcore/test/testRooGradMinimizerFcn.cxx b/roofit/roofitcore/test/testRooGradMinimizerFcn.cxx
index 3d1a8181612c8..d1ee53aec994e 100644
--- a/roofit/roofitcore/test/testRooGradMinimizerFcn.cxx
+++ b/roofit/roofitcore/test/testRooGradMinimizerFcn.cxx
@@ -29,69 +29,73 @@
 #include <RooMsgService.h>
 #include <RooGlobalFunc.h> // RooFit::ERROR
 
-TEST(GradMinimizer, Gaussian1D)
+class GradMinimizerParSeed : public testing::TestWithParam<unsigned long> {};
+
+TEST_P(GradMinimizerParSeed, Gaussian1D)
 {
    RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
 
    // produce the same random stuff every time
-   RooRandom::randomGenerator()->SetSeed(1);
+   RooRandom::randomGenerator()->SetSeed(GetParam());
 
-   for (int i = 0; i < 10; ++i) {
-      RooWorkspace w = RooWorkspace();
-
-      std::unique_ptr<RooAbsReal> nll;
-      std::unique_ptr<RooArgSet> values;
-      RooAbsPdf *pdf;
-      RooDataSet *data;
-      std::tie(nll, pdf, data, values) = generate_1D_gaussian_pdf_nll(w, 10000);
-      // when c++17 support arrives, change to this:
-      //  auto [nll, values] = generate_1D_gaussian_pdf_nll(w, 10000);
-      RooRealVar *mu = w.var("mu");
-
-      RooArgSet *savedValues = dynamic_cast<RooArgSet *>(values->snapshot());
-      if (savedValues == nullptr) {
-         throw std::runtime_error("params->snapshot() cannot be casted to RooArgSet!");
-      }
+   RooWorkspace w = RooWorkspace();
 
-      // --------
+   std::unique_ptr<RooAbsReal> nll;
+   std::unique_ptr<RooArgSet> values;
+   RooAbsPdf *pdf;
+   RooDataSet *data;
+   std::tie(nll, pdf, data, values) = generate_1D_gaussian_pdf_nll(w, 10000);
+   // when c++17 support arrives, change to this:
+   //  auto [nll, values] = generate_1D_gaussian_pdf_nll(w, 10000);
+   RooRealVar *mu = w.var("mu");
 
-      RooMinimizer m0(*nll);
-      m0.setMinimizerType("Minuit2");
+   RooArgSet *savedValues = dynamic_cast<RooArgSet *>(values->snapshot());
+   if (savedValues == nullptr) {
+      throw std::runtime_error("params->snapshot() cannot be casted to RooArgSet!");
+   }
+
+   // --------
 
-      m0.setStrategy(0);
-      m0.setPrintLevel(-1);
+   RooMinimizer m0(*nll);
+   m0.setMinimizerType("Minuit2");
+
+   m0.setStrategy(0);
+   m0.setPrintLevel(-1);
 
-      m0.migrad();
+   m0.migrad();
 
-      RooFitResult *m0result = m0.lastMinuitFit();
-      double minNll0 = m0result->minNll();
-      double edm0 = m0result->edm();
-      double mu0 = mu->getVal();
-      double muerr0 = mu->getError();
+   RooFitResult *m0result = m0.lastMinuitFit();
+   double minNll0 = m0result->minNll();
+   double edm0 = m0result->edm();
+   double mu0 = mu->getVal();
+   double muerr0 = mu->getError();
 
-      *values = *savedValues;
+   *values = *savedValues;
 
-      std::unique_ptr<RooMinimizer> m1 = RooMinimizer::create(*nll, RooMinimizer::FcnMode::gradient);
-      m1->setMinimizerType("Minuit2");
+   std::unique_ptr<RooMinimizer> m1 = RooMinimizer::create(*nll, RooMinimizer::FcnMode::gradient);
+   m1->setMinimizerType("Minuit2");
 
-      m1->setStrategy(0);
-      m1->setPrintLevel(-1);
+   m1->setStrategy(0);
+   m1->setPrintLevel(-1);
 
-      m1->migrad();
+   m1->migrad();
 
-      RooFitResult *m1result = m1->lastMinuitFit();
-      double minNll1 = m1result->minNll();
-      double edm1 = m1result->edm();
-      double mu1 = mu->getVal();
-      double muerr1 = mu->getError();
+   RooFitResult *m1result = m1->lastMinuitFit();
+   double minNll1 = m1result->minNll();
+   double edm1 = m1result->edm();
+   double mu1 = mu->getVal();
+   double muerr1 = mu->getError();
 
-      EXPECT_EQ(minNll0, minNll1);
-      EXPECT_EQ(mu0, mu1);
-      EXPECT_EQ(muerr0, muerr1);
-      EXPECT_EQ(edm0, edm1);
-   }
+   EXPECT_EQ(minNll0, minNll1);
+   EXPECT_EQ(mu0, mu1);
+   EXPECT_EQ(muerr0, muerr1);
+   EXPECT_EQ(edm0, edm1);
 }
 
+INSTANTIATE_TEST_SUITE_P(Seeds,
+                         GradMinimizerParSeed,
+                         testing::Range<unsigned long>(1, 11));
+
 TEST(GradMinimizerDebugging, DISABLED_Gaussian1DNominal)
 {
    // produce the same random stuff every time
@@ -182,8 +186,8 @@ TEST(GradMinimizer, GaussianND)
    RooFitResult *m0result = m0.lastMinuitFit();
    double minNll0 = m0result->minNll();
    double edm0 = m0result->edm();
-   double mean0[N];
-   double std0[N];
+   std::vector<double> mean0(N);
+   std::vector<double> std0(N);
    for (unsigned ix = 0; ix < N; ++ix) {
       {
          std::ostringstream os;
@@ -214,8 +218,8 @@ TEST(GradMinimizer, GaussianND)
    RooFitResult *m1result = m1->lastMinuitFit();
    double minNll1 = m1result->minNll();
    double edm1 = m1result->edm();
-   double mean1[N];
-   double std1[N];
+   std::vector<double> mean1(N);
+   std::vector<double> std1(N);
    for (unsigned ix = 0; ix < N; ++ix) {
       {
          std::ostringstream os;
@@ -279,8 +283,8 @@ TEST(GradMinimizerReverse, GaussianND)
    RooFitResult *m0result = m0->lastMinuitFit();
    double minNll0 = m0result->minNll();
    double edm0 = m0result->edm();
-   double mean0[N];
-   double std0[N];
+   std::vector<double> mean0(N);
+   std::vector<double> std0(N);
    for (unsigned ix = 0; ix < N; ++ix) {
       {
          std::ostringstream os;
@@ -311,8 +315,8 @@ TEST(GradMinimizerReverse, GaussianND)
    RooFitResult *m1result = m1.lastMinuitFit();
    double minNll1 = m1result->minNll();
    double edm1 = m1result->edm();
-   double mean1[N];
-   double std1[N];
+   std::vector<double> mean1(N);
+   std::vector<double> std1(N);
    for (unsigned ix = 0; ix < N; ++ix) {
       {
          std::ostringstream os;
diff --git a/roofit/roofitcore/test/test_lib.h b/roofit/roofitcore/test/test_lib.h
index dae7ef3b50dd7..df6ead82c263d 100644
--- a/roofit/roofitcore/test/test_lib.h
+++ b/roofit/roofitcore/test/test_lib.h
@@ -21,6 +21,7 @@
 
 #include <sstream>
 #include <memory>  // make_unique
+#include <vector>
 
 RooAbsPdf * generate_1D_gaussian_pdf(RooWorkspace &w)
 {
@@ -61,7 +62,7 @@ generate_ND_gaussian_pdf_nll(RooWorkspace &w, unsigned int n, unsigned long N_ev
   RooArgSet obs_set;
 
   // create gaussian parameters
-  double mean[n], sigma[n];
+  std::vector<double> mean(n), sigma(n);
   for (unsigned ix = 0; ix < n; ++ix) {
     mean[ix] = RooRandom::randomGenerator()->Gaus(0, 2);
     sigma[ix] = 0.1 + abs(RooRandom::randomGenerator()->Gaus(0, 2));