diff --git a/test/ExampleProblems/BellCurveLeastSquaresProblem.cs b/test/ExampleProblems/CustomProblems/BellCurveProblem.cs
similarity index 96%
rename from test/ExampleProblems/BellCurveLeastSquaresProblem.cs
rename to test/ExampleProblems/CustomProblems/BellCurveProblem.cs
index 2dc2463..510ace3 100644
--- a/test/ExampleProblems/BellCurveLeastSquaresProblem.cs
+++ b/test/ExampleProblems/CustomProblems/BellCurveProblem.cs
@@ -1,6 +1,6 @@
-namespace ExampleProblems;
+namespace ExampleProblems.CustomProblems;
 
-public class BellCurveLeastSquaresProblem : ConfigurableLeastSquaresProblem
+public class BellCurveProblem : ConfigurableLeastSquaresProblem
 {
     private static double Bell(double x, double location, double variance)
     {
diff --git a/test/ExampleProblems/CubicPolynomialLeastSquaresProblem.cs b/test/ExampleProblems/CustomProblems/CubicPolynomialProblem.cs
similarity index 93%
rename from test/ExampleProblems/CubicPolynomialLeastSquaresProblem.cs
rename to test/ExampleProblems/CustomProblems/CubicPolynomialProblem.cs
index 5a58c40..47834fa 100644
--- a/test/ExampleProblems/CubicPolynomialLeastSquaresProblem.cs
+++ b/test/ExampleProblems/CustomProblems/CubicPolynomialProblem.cs
@@ -1,6 +1,6 @@
-namespace ExampleProblems;
+namespace ExampleProblems.CustomProblems;
 
-public class CubicPolynomialLeastSquaresProblem : ConfigurableLeastSquaresProblem
+public class CubicPolynomialProblem : ConfigurableLeastSquaresProblem
 {
     protected override Func<double, IReadOnlyList<double>, double> Model { get; } =
         (x, c) => c[0] + c[1] * x + c[2] * x * x + c[3] * x * x * x;
diff --git a/test/ExampleProblems/ExponentialDecayLeastSquaresProblem.cs b/test/ExampleProblems/CustomProblems/ExponentialDecayProblem.cs
similarity index 94%
rename from test/ExampleProblems/ExponentialDecayLeastSquaresProblem.cs
rename to test/ExampleProblems/CustomProblems/ExponentialDecayProblem.cs
index 35312ea..f30b180 100644
--- a/test/ExampleProblems/ExponentialDecayLeastSquaresProblem.cs
+++ b/test/ExampleProblems/CustomProblems/ExponentialDecayProblem.cs
@@ -1,6 +1,6 @@
-namespace ExampleProblems;
+namespace ExampleProblems.CustomProblems;
 
-public class ExponentialDecayLeastSquaresProblem : ConfigurableLeastSquaresProblem
+public class ExponentialDecayProblem : ConfigurableLeastSquaresProblem
 {
     protected override Func<double, IReadOnlyList<double>, double> Model { get; } =
         (x, p) => p[0] * Math.Exp(-p[1] * x) + p[2];
diff --git a/test/ExampleProblems/NumericalPendulumLeastSquaresProblem.cs b/test/ExampleProblems/CustomProblems/NumericalPendulumProblem.cs
similarity index 78%
rename from test/ExampleProblems/NumericalPendulumLeastSquaresProblem.cs
rename to test/ExampleProblems/CustomProblems/NumericalPendulumProblem.cs
index aca027e..490ed3d 100644
--- a/test/ExampleProblems/NumericalPendulumLeastSquaresProblem.cs
+++ b/test/ExampleProblems/CustomProblems/NumericalPendulumProblem.cs
@@ -1,17 +1,17 @@
 using minuit2.net;
 using minuit2.net.CostFunctions;
 
-namespace ExampleProblems;
+namespace ExampleProblems.CustomProblems;
 
-public class NumericalPendulumLeastSquaresProblem : IConfiguredProblem
+public class NumericalPendulumProblem : IConfiguredProblem
 {
-    public NumericalPendulumLeastSquaresProblem()
+    public NumericalPendulumProblem()
     {
         OptimumParameterValues = [1];
         ParameterConfigurations = [ParameterConfiguration.Variable("pendulumLength", 2)];
 
         var model = NumericalPendulumModelFor(startAngle: 1.5, startAngleVelocity: 0.0);
-        var xValues = LinearlySpacedValues(0, 3, 0.001);
+        var xValues = Values.LinearlySpacedBetween(0, 3, 0.001);
         var yValues = model(xValues, OptimumParameterValues.ToArray());
         Cost = CostFunction.LeastSquares(xValues, yValues, ["pendulumLength"], model);
     }
@@ -40,10 +40,4 @@ private static Func<IReadOnlyList<double>, IReadOnlyList<double>, IReadOnlyList<
             return angle;
         };
     }
-
-    private static double[] LinearlySpacedValues(double start, double end, double step)
-    {
-        var number = (int)((end - start) / step) + 1;
-        return Enumerable.Range(0, number).Select(i => start + i * step).ToArray();
-    }
 }
\ No newline at end of file
diff --git a/test/ExampleProblems/QuadraticPolynomialLeastSquaresProblem.cs b/test/ExampleProblems/CustomProblems/QuadraticPolynomialProblem.cs
similarity index 93%
rename from test/ExampleProblems/QuadraticPolynomialLeastSquaresProblem.cs
rename to test/ExampleProblems/CustomProblems/QuadraticPolynomialProblem.cs
index 423695d..c598c14 100644
--- a/test/ExampleProblems/QuadraticPolynomialLeastSquaresProblem.cs
+++ b/test/ExampleProblems/CustomProblems/QuadraticPolynomialProblem.cs
@@ -1,6 +1,6 @@
-namespace ExampleProblems;
+namespace ExampleProblems.CustomProblems;
 
-public class QuadraticPolynomialLeastSquaresProblem : ConfigurableLeastSquaresProblem
+public class QuadraticPolynomialProblem : ConfigurableLeastSquaresProblem
 {
     protected override Func<double, IReadOnlyList<double>, double> Model { get; } =
         (x, c) => c[0] + c[1] * x + c[2] * x * x;
diff --git a/test/ExampleProblems/CustomProblems/SurfaceBiosensorBindingKineticsProblem.cs b/test/ExampleProblems/CustomProblems/SurfaceBiosensorBindingKineticsProblem.cs
new file mode 100644
index 0000000..9d94860
--- /dev/null
+++ b/test/ExampleProblems/CustomProblems/SurfaceBiosensorBindingKineticsProblem.cs
@@ -0,0 +1,245 @@
+using minuit2.net;
+using minuit2.net.CostFunctions;
+using static ExampleProblems.DerivativeConfiguration;
+using static minuit2.net.ParameterConfiguration;
+
+namespace ExampleProblems.CustomProblems;
+
+public class SurfaceBiosensorBindingKineticsProblem : IConfiguredProblem
+{
+   // Simplified assumptions:
+   // - ligand is immobilized on surface with surface density remaining constant
+   // - analyte association kinetics (starting at t = 0) is directly followed by analyte dissociation kinetics
+   // - analyte concentration is 0 during the dissociation measurement (no rebinding)
+   // - analyte concentration remains constant during the association measurement
+   // - analyte mobility is no restricting factor for the association kinetics (no mass transport limitations)
+   
+   private readonly DerivativeConfiguration _modelDerivativeConfiguration;
+   private readonly IReadOnlyList<double> _x;
+   private readonly IReadOnlyList<double> _y;
+   private readonly IReadOnlyList<string> _parameters;
+
+   public SurfaceBiosensorBindingKineticsProblem(
+      DerivativeConfiguration modelDerivativeConfiguration = WithoutDerivatives,
+      double analyteConcentrationValueInNanoMolar = 10,
+      string? localParameterSuffix = null,
+      int randomSeed = 0)
+   {
+      _modelDerivativeConfiguration = modelDerivativeConfiguration;
+      var random = new Random(randomSeed);
+
+      var amplitude = Local("amp");
+      const string associationRate = "ka";
+      var analyteConcentration = Local("c");
+      const string dissociationRate = "kd";
+      var dissociationStart = Local("td");
+      
+      _parameters = [amplitude, associationRate, analyteConcentration, dissociationRate, dissociationStart];
+      var parameterValues = ParameterValuesFor(analyteConcentrationValueInNanoMolar);
+
+      _x = Values.LinearlySpacedBetween(0, 1000, 1);
+      _y = _x.Select(x => Model(x, parameterValues) + random.NextNormal(0, 0.001)).ToArray();
+
+      OptimumParameterValues = parameterValues;
+      ParameterConfigurations =
+      [
+         Variable(amplitude, parameterValues[0] * random.NextNormal(1, 0.1)),
+         Variable(associationRate, parameterValues[1] * random.NextNormal(1, 0.1), lowerLimit: 0),
+         Fixed(analyteConcentration, parameterValues[2]),
+         Variable(dissociationRate, parameterValues[3] * random.NextNormal(1, 0.1), lowerLimit: 0),
+         Variable(dissociationStart, parameterValues[4] + random.NextNormal(0, 3))
+      ];
+      return;
+
+      string Local(string parameter) => localParameterSuffix == null ? parameter : parameter + localParameterSuffix;
+   }
+   
+   public static IConfiguredProblem Global(
+      IEnumerable<double> analyteConcentrationValuesInNanoMolar,
+      DerivativeConfiguration modelDerivativeConfiguration = WithoutDerivatives,
+      int randomSeed = 0)
+   {
+      var problems = analyteConcentrationValuesInNanoMolar.Select((c, i) =>
+         new SurfaceBiosensorBindingKineticsProblem(modelDerivativeConfiguration, c, i.ToString(), i + randomSeed));
+      
+      var costs = new List<ICostFunction>();
+      var parameterConfigs = new List<ParameterConfiguration>();
+      var optimumValues = new List<double>();
+      foreach (var problem in problems)
+      {
+         costs.Add(problem.Cost);
+         foreach (var (config, value) in problem.ParameterConfigurations.Zip(problem.OptimumParameterValues))
+         {
+            if (parameterConfigs.Any(c => c.Name == config.Name)) continue;
+            parameterConfigs.Add(config);
+            optimumValues.Add(value);
+         }
+      }
+
+      return new ConfiguredProblem(CostFunction.Sum(costs.ToArray()), optimumValues, parameterConfigs);
+   }
+
+   private static IReadOnlyList<double> ParameterValuesFor(double analyteConcentrationInNanoMolar)
+   {
+      const double ka = 1e-3; // (s * nM)-1
+      const double kd = 1e-2; // s-1
+      var kaObs = ka * analyteConcentrationInNanoMolar + kd;
+      var amp = ka * analyteConcentrationInNanoMolar / kaObs;
+      var td = Math.Max(5 / kaObs, 100);
+      return [amp, ka, analyteConcentrationInNanoMolar, kd, td];
+   }
+
+   private static Func<double, IReadOnlyList<double>, double> Model { get; } = (x, p) =>
+   {
+      var (amp, ka, c, kd, td) = (p[0], p[1], p[2], p[3], p[4]);
+      var kaObs = ka * c + kd;
+
+      return td > x
+         ? amp * (1 - Math.Exp(-kaObs * x))
+         : amp * (1 - Math.Exp(-kaObs * td)) * Math.Exp(-kd * (x - td));
+   };
+
+   private static Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelGradient { get; } = (x, p) =>
+   {
+      var (amp, ka, c, kd, td) = (p[0], p[1], p[2], p[3], p[4]);
+      var kaObs = ka * c + kd;
+
+      var g0 = td > x ? 1 - Math.Exp(-kaObs * x) : (1 - Math.Exp(-kaObs * td)) * Math.Exp(-kd * (x - td));
+      var g1 = td > x
+         ? amp * c * x * Math.Exp(-kaObs * x)
+         : amp * c * td * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var g2 = td > x
+         ? amp * ka * x * Math.Exp(-kaObs * x)
+         : amp * ka * td * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var g3 = td > x
+         ? amp * x * Math.Exp(-kaObs * x)
+         : amp * td * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) +
+           amp * (1 - Math.Exp(-kaObs * td)) * (td - x) * Math.Exp(-kd * (x - td));
+      var g4 = td > x
+         ? 0
+         : amp * kd * (1 - Math.Exp(-kaObs * td)) * Math.Exp(-kd * (x - td)) -
+           amp * -kaObs * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+
+      return [g0, g1, g2, g3, g4];
+   };
+
+   private static Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessian { get; } = (x, p) =>
+   {
+      var (amp, ka, c, kd, td) = (p[0], p[1], p[2], p[3], p[4]);
+      var kaObs = ka * c + kd;
+
+      const double h00 = 0;
+      var h01 = td > x
+         ? c * x * Math.Exp(-kaObs * x)
+         : c * td * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var h02 = td > x
+         ? ka * x * Math.Exp(-kaObs * x)
+         : ka * td * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var h03 = td > x
+         ? x * Math.Exp(-kaObs * x)
+         : td * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) +
+           (1 - Math.Exp(-kaObs * td)) * (td - x) * Math.Exp(-kd * (x - td));
+      var h04 = td > x
+         ? 0
+         : kd * (1 - Math.Exp(-kaObs * td)) * Math.Exp(-kd * (x - td)) -
+           -kaObs * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var h11 = td > x
+         ? -amp * Math.Pow(c, 2) * Math.Pow(x, 2) * Math.Exp(-kaObs * x)
+         : -amp * Math.Pow(c, 2) * Math.Pow(td, 2) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var h12 = td > x
+         ? -amp * c * ka * Math.Pow(x, 2) * Math.Exp(-kaObs * x) + amp * x * Math.Exp(-kaObs * x)
+         : -amp * c * ka * Math.Pow(td, 2) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) +
+           amp * td * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var h13 = td > x
+         ? -amp * c * Math.Pow(x, 2) * Math.Exp(-kaObs * x)
+         : -amp * c * Math.Pow(td, 2) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) +
+           amp * c * td * (td - x) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var h14 = td > x
+         ? 0
+         : amp * c * kd * td * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) +
+           amp * c * -kaObs * td * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) +
+           amp * c * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var h22 = td > x
+         ? -amp * Math.Pow(ka, 2) * Math.Pow(x, 2) * Math.Exp(-kaObs * x)
+         : -amp * Math.Pow(ka, 2) * Math.Pow(td, 2) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var h23 = td > x
+         ? -amp * ka * Math.Pow(x, 2) * Math.Exp(-kaObs * x)
+         : -amp * ka * Math.Pow(td, 2) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) +
+           amp * ka * td * (td - x) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var h24 = td > x
+         ? 0
+         : amp * ka * kd * td * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) +
+           amp * ka * -kaObs * td * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) +
+           amp * ka * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var h33 = td > x
+         ? -amp * Math.Pow(x, 2) * Math.Exp(-kaObs * x)
+         : -amp * Math.Pow(td, 2) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) +
+           2 * amp * td * (td - x) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) + amp *
+           (1 - Math.Exp(-kaObs * td)) * Math.Pow(td - x, 2) * Math.Exp(-kd * (x - td));
+      var h34 = td > x
+         ? 0
+         : amp * kd * td * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) +
+           amp * kd * (1 - Math.Exp(-kaObs * td)) * (td - x) * Math.Exp(-kd * (x - td)) +
+           amp * -kaObs * td * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) +
+           amp * (1 - Math.Exp(-kaObs * td)) * Math.Exp(-kd * (x - td)) -
+           amp * (td - x) * -kaObs * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) +
+           amp * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var h44 = td > x
+         ? 0
+         : amp * Math.Pow(kd, 2) * (1 - Math.Exp(-kaObs * td)) * Math.Exp(-kd * (x - td)) -
+           2 * amp * kd * -kaObs * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) -
+           amp * Math.Pow(-kaObs, 2) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+
+      return
+      [
+         h00, h01, h02, h03, h04,
+         h01, h11, h12, h13, h14,
+         h02, h12, h22, h23, h24,
+         h03, h13, h23, h33, h34,
+         h04, h14, h24, h34, h44,
+      ];
+   };
+
+   private static Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessianDiagonal { get; } = (x, p) =>
+   {
+      var (amp, ka, c, kd, td) = (p[0], p[1], p[2], p[3], p[4]);
+      var kaObs = ka * c + kd;
+
+      const double h00 = 0;
+      var h11 = td > x
+         ? -amp * Math.Pow(c, 2) * Math.Pow(x, 2) * Math.Exp(-kaObs * x)
+         : -amp * Math.Pow(c, 2) * Math.Pow(td, 2) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var h22 = td > x
+         ? -amp * Math.Pow(ka, 2) * Math.Pow(x, 2) * Math.Exp(-kaObs * x)
+         : -amp * Math.Pow(ka, 2) * Math.Pow(td, 2) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+      var h33 = td > x
+         ? -amp * Math.Pow(x, 2) * Math.Exp(-kaObs * x)
+         : -amp * Math.Pow(td, 2) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) +
+           2 * amp * td * (td - x) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) + amp *
+           (1 - Math.Exp(-kaObs * td)) * Math.Pow(td - x, 2) * Math.Exp(-kd * (x - td));
+      var h44 = td > x
+         ? 0
+         : amp * Math.Pow(kd, 2) * (1 - Math.Exp(-kaObs * td)) * Math.Exp(-kd * (x - td)) -
+           2 * amp * kd * -kaObs * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td) -
+           amp * Math.Pow(-kaObs, 2) * Math.Exp(-kd * (x - td)) * Math.Exp(-kaObs * td);
+
+      return [h00, h11, h22, h33, h44];
+   };
+
+   public ICostFunction Cost => _modelDerivativeConfiguration switch
+   {
+      WithoutDerivatives => 
+         CostFunction.LeastSquares(_x, _y, _parameters, Model),
+      WithGradient => 
+         CostFunction.LeastSquares(_x, _y, _parameters, Model, ModelGradient),
+      WithGradientAndHessian => 
+         CostFunction.LeastSquares(_x, _y, _parameters, Model, ModelGradient, ModelHessian),
+      WithGradientHessianAndHessianDiagonal => 
+         CostFunction.LeastSquares(_x, _y, _parameters, Model, ModelGradient, ModelHessian, ModelHessianDiagonal),
+      _ => throw new ArgumentOutOfRangeException(nameof(_modelDerivativeConfiguration), _modelDerivativeConfiguration, null)
+   };
+
+   public IReadOnlyCollection<double> OptimumParameterValues { get; }
+
+   public IReadOnlyCollection<ParameterConfiguration> ParameterConfigurations { get; }
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/DerivativeConfiguration.cs b/test/ExampleProblems/DerivativeConfiguration.cs
new file mode 100644
index 0000000..363368f
--- /dev/null
+++ b/test/ExampleProblems/DerivativeConfiguration.cs
@@ -0,0 +1,9 @@
+namespace ExampleProblems;
+
+public enum DerivativeConfiguration
+{
+    WithoutDerivatives,
+    WithGradient,
+    WithGradientAndHessian,
+    WithGradientHessianAndHessianDiagonal
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/DirectCostFunction.cs b/test/ExampleProblems/DirectCostFunction.cs
new file mode 100644
index 0000000..ea30390
--- /dev/null
+++ b/test/ExampleProblems/DirectCostFunction.cs
@@ -0,0 +1,25 @@
+using minuit2.net;
+using minuit2.net.CostFunctions;
+
+namespace ExampleProblems;
+
+internal class DirectCostFunction(
+    IReadOnlyList<string> parameters,
+    Func<IReadOnlyList<double>, double> value, 
+    Func<IReadOnlyList<double>, IReadOnlyList<double>>? gradient = null, 
+    Func<IReadOnlyList<double>, IReadOnlyList<double>>? hessian = null, 
+    Func<IReadOnlyList<double>, IReadOnlyList<double>>? hessianDiagonal = null, 
+    double errorDefinition = 1) 
+    : ICostFunction
+{
+    public IReadOnlyList<string> Parameters { get; } = parameters;
+    public bool HasGradient { get; } = gradient != null;
+    public bool HasHessian { get; } = hessian != null;
+    public bool HasHessianDiagonal { get; } = hessianDiagonal != null;
+    public double ErrorDefinition { get; } = errorDefinition;
+    public double ValueFor(IReadOnlyList<double> parameterValues) => value(parameterValues);
+    public IReadOnlyList<double> GradientFor(IReadOnlyList<double> parameterValues) => gradient!(parameterValues);
+    public IReadOnlyList<double> HessianFor(IReadOnlyList<double> parameterValues) => hessian!(parameterValues);
+    public IReadOnlyList<double> HessianDiagonalFor(IReadOnlyList<double> parameterValues) => hessianDiagonal!(parameterValues);
+    public ICostFunction WithErrorDefinitionRecalculatedBasedOnValid(IMinimizationResult result) => this;
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/MinuitTutorialProblems/FletcherPowellProblem.cs b/test/ExampleProblems/MinuitTutorialProblems/FletcherPowellProblem.cs
new file mode 100644
index 0000000..9e9a423
--- /dev/null
+++ b/test/ExampleProblems/MinuitTutorialProblems/FletcherPowellProblem.cs
@@ -0,0 +1,155 @@
+namespace ExampleProblems.MinuitTutorialProblems;
+
+public class FletcherPowellProblem(DerivativeConfiguration derivativeConfiguration)
+    : MinuitTutorialProblem(
+        Parameters,
+        InitialValues,
+        OptimumValues,
+        Function,
+        Gradient,
+        Hessian,
+        HessianDiagonal,
+        derivativeConfiguration)
+{
+    // see section 7.4 in https://seal.web.cern.ch/seal/documents/minuit/mntutorial.pdf
+    // NOTE: start and optimum parameter values are swapped in the source
+
+    private static readonly IReadOnlyList<string> Parameters = ["x", "y", "z"];
+    private static readonly IReadOnlyList<double> InitialValues = [1, 0, 0];
+    private static readonly IReadOnlyList<double> OptimumValues = [-1, 0, 0];
+
+    private static readonly Func<IReadOnlyList<double>, double> Function = p =>
+    {
+        var (x, y, z) = (p[0], p[1], p[2]);
+        return x < 0
+            ? Math.Pow(z, 2) + 100 * Math.Pow(z - 5 * Math.Atan(y / x) / Math.PI, 2) +
+              100 * Math.Pow(Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) - 1, 2)
+            : Math.Pow(z, 2) + 100 * Math.Pow(z - 1d / 2 * (10 * Math.Atan(y / x) + 10 * Math.PI) / Math.PI, 2) +
+              100 * Math.Pow(Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) - 1, 2);
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> Gradient = p =>
+    {
+        var (x, y, z) = (p[0], p[1], p[2]);
+        var g0 = x < 0
+            ? 200 *
+              (Math.Pow(Math.PI, 2) * x * (Math.Pow(x, 2) + Math.Pow(y, 2)) *
+               (Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) - 1) +
+               5 * y * Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) * (Math.PI * z - 5 * Math.Atan(y / x))) /
+              (Math.Pow(Math.PI, 2) * Math.Pow(Math.Pow(x, 2) + Math.Pow(y, 2), 3d / 2))
+            : 200 *
+              (Math.Pow(Math.PI, 2) * x * (Math.Pow(x, 2) + Math.Pow(y, 2)) *
+               (Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) - 1) +
+               5 * y * Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) *
+               (Math.PI * z - 5 * Math.Atan(y / x) - 5 * Math.PI)) /
+              (Math.Pow(Math.PI, 2) * Math.Pow(Math.Pow(x, 2) + Math.Pow(y, 2), 3d / 2));
+        var g1 = x < 0
+            ? 200 * (-5 * x * Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) * (Math.PI * z - 5 * Math.Atan(y / x)) +
+                     Math.Pow(Math.PI, 2) * y * (Math.Pow(x, 2) + Math.Pow(y, 2)) *
+                     (Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) - 1)) /
+              (Math.Pow(Math.PI, 2) * Math.Pow(Math.Pow(x, 2) + Math.Pow(y, 2), 3d / 2))
+            : 200 *
+              (-5 * x * Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) *
+               (Math.PI * z - 5 * Math.Atan(y / x) - 5 * Math.PI) +
+               Math.Pow(Math.PI, 2) * y * (Math.Pow(x, 2) + Math.Pow(y, 2)) *
+               (Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) - 1)) /
+              (Math.Pow(Math.PI, 2) * Math.Pow(Math.Pow(x, 2) + Math.Pow(y, 2), 3d / 2));
+        var g2 = x < 0
+            ? 202 * z - 1000 * Math.Atan(y / x) / Math.PI
+            : 202 * z - 1000 * Math.Atan(y / x) / Math.PI - 1000;
+        return [g0, g1, g2];
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> Hessian = p =>
+    {
+        var (x, y, z) = (p[0], p[1], p[2]);
+        var h00 = x < 0
+            ? 200 * (Math.Pow(Math.PI, 2) * Math.Pow(x, 4) +
+                     2 * Math.Pow(Math.PI, 2) * Math.Pow(x, 2) * Math.Pow(y, 2) - 10 * Math.PI * x * y * z +
+                     50 * x * y * Math.Atan(y / x) + Math.Pow(Math.PI, 2) * Math.Pow(y, 4) -
+                     Math.Pow(Math.PI, 2) * Math.Pow(y, 2) * Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) +
+                     25 * Math.Pow(y, 2)) /
+              (Math.Pow(Math.PI, 2) * (Math.Pow(x, 4) + 2 * Math.Pow(x, 2) * Math.Pow(y, 2) + Math.Pow(y, 4)))
+            : 200 * (Math.Pow(Math.PI, 2) * Math.Pow(x, 4) +
+                     2 * Math.Pow(Math.PI, 2) * Math.Pow(x, 2) * Math.Pow(y, 2) - 10 * Math.PI * x * y * z +
+                     50 * x * y * Math.Atan(y / x) + 50 * Math.PI * x * y + Math.Pow(Math.PI, 2) * Math.Pow(y, 4) -
+                     Math.Pow(Math.PI, 2) * Math.Pow(y, 2) * Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) +
+                     25 * Math.Pow(y, 2)) / (Math.Pow(Math.PI, 2) *
+                                             (Math.Pow(x, 4) + 2 * Math.Pow(x, 2) * Math.Pow(y, 2) + Math.Pow(y, 4)));
+        var h01 = x < 0
+            ? -5000 * x * y / (Math.Pow(Math.PI, 2) * Math.Pow(Math.Pow(x, 2) + Math.Pow(y, 2), 2)) +
+              200 * x * y / Math.Pow(Math.Pow(x, 2) + Math.Pow(y, 2), 3d / 2) -
+              2000 * Math.Pow(y, 2) * z / (Math.PI * Math.Pow(Math.Pow(x, 2) + Math.Pow(y, 2), 2)) +
+              10000 * Math.Pow(y, 2) * Math.Atan(y / x) /
+              (Math.Pow(Math.PI, 2) * Math.Pow(Math.Pow(x, 2) + Math.Pow(y, 2), 2)) +
+              1000 * z / (Math.PI * Math.Pow(x, 2) + Math.PI * Math.Pow(y, 2)) -
+              5000 * Math.Atan(y / x) / (Math.Pow(Math.PI, 2) * Math.Pow(x, 2) + Math.Pow(Math.PI, 2) * Math.Pow(y, 2))
+            : 200 * (Math.Pow(Math.PI, 2) * x * y * Math.Pow(Math.Pow(x, 2) + Math.Pow(y, 2), 9d / 2) -
+                     Math.Pow(Math.PI, 2) * x * y * Math.Pow(Math.Pow(x, 2) + Math.Pow(y, 2), 4) *
+                     (Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) - 1) +
+                     5 * y * (-5 * x + 2 * y * (-Math.PI * z + 5 * Math.Atan(y / x) + 5 * Math.PI)) *
+                     Math.Pow(Math.Pow(x, 2) + Math.Pow(y, 2), 7d / 2) +
+                     5 * Math.Pow(Math.Pow(x, 2) + Math.Pow(y, 2), 9d / 2) *
+                     (Math.PI * z - 5 * Math.Atan(y / x) - 5 * Math.PI)) /
+              (Math.Pow(Math.PI, 2) * Math.Pow(Math.Pow(x, 2) + Math.Pow(y, 2), 11d / 2));
+        var h02 = 1000 * y / (Math.PI * (Math.Pow(x, 2) + Math.Pow(y, 2)));
+        var h11 = x < 0
+            ? 200 * (Math.Pow(Math.PI, 2) * Math.Pow(x, 4) +
+                     2 * Math.Pow(Math.PI, 2) * Math.Pow(x, 2) * Math.Pow(y, 2) -
+                     Math.Pow(Math.PI, 2) * Math.Pow(x, 2) * Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) +
+                     25 * Math.Pow(x, 2) + 10 * Math.PI * x * y * z - 50 * x * y * Math.Atan(y / x) +
+                     Math.Pow(Math.PI, 2) * Math.Pow(y, 4)) /
+              (Math.Pow(Math.PI, 2) * (Math.Pow(x, 4) + 2 * Math.Pow(x, 2) * Math.Pow(y, 2) + Math.Pow(y, 4)))
+            : 200 * (Math.Pow(Math.PI, 2) * Math.Pow(x, 4) +
+                     2 * Math.Pow(Math.PI, 2) * Math.Pow(x, 2) * Math.Pow(y, 2) -
+                     Math.Pow(Math.PI, 2) * Math.Pow(x, 2) * Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) +
+                     25 * Math.Pow(x, 2) + 10 * Math.PI * x * y * z - 50 * x * y * Math.Atan(y / x) -
+                     50 * Math.PI * x * y +
+                     Math.Pow(Math.PI, 2) * Math.Pow(y, 4)) / (Math.Pow(Math.PI, 2) *
+                                                               (Math.Pow(x, 4) + 2 * Math.Pow(x, 2) * Math.Pow(y, 2) +
+                                                                Math.Pow(y, 4)));
+        var h12 = -1000 * x / (Math.PI * (Math.Pow(x, 2) + Math.Pow(y, 2)));
+        const double h22 = 202;
+        return
+        [
+            h00, h01, h02,
+            h01, h11, h12,
+            h02, h12, h22
+        ];
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> HessianDiagonal = p =>
+    {
+        var (x, y, z) = (p[0], p[1], p[2]);
+        var h00 = x < 0
+            ? 200 * (Math.Pow(Math.PI, 2) * Math.Pow(x, 4) +
+                     2 * Math.Pow(Math.PI, 2) * Math.Pow(x, 2) * Math.Pow(y, 2) - 10 * Math.PI * x * y * z +
+                     50 * x * y * Math.Atan(y / x) + Math.Pow(Math.PI, 2) * Math.Pow(y, 4) -
+                     Math.Pow(Math.PI, 2) * Math.Pow(y, 2) * Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) +
+                     25 * Math.Pow(y, 2)) /
+              (Math.Pow(Math.PI, 2) * (Math.Pow(x, 4) + 2 * Math.Pow(x, 2) * Math.Pow(y, 2) + Math.Pow(y, 4)))
+            : 200 * (Math.Pow(Math.PI, 2) * Math.Pow(x, 4) +
+                     2 * Math.Pow(Math.PI, 2) * Math.Pow(x, 2) * Math.Pow(y, 2) - 10 * Math.PI * x * y * z +
+                     50 * x * y * Math.Atan(y / x) + 50 * Math.PI * x * y + Math.Pow(Math.PI, 2) * Math.Pow(y, 4) -
+                     Math.Pow(Math.PI, 2) * Math.Pow(y, 2) * Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) +
+                     25 * Math.Pow(y, 2)) / (Math.Pow(Math.PI, 2) *
+                                             (Math.Pow(x, 4) + 2 * Math.Pow(x, 2) * Math.Pow(y, 2) + Math.Pow(y, 4)));
+        var h11 = x < 0
+            ? 200 * (Math.Pow(Math.PI, 2) * Math.Pow(x, 4) +
+                     2 * Math.Pow(Math.PI, 2) * Math.Pow(x, 2) * Math.Pow(y, 2) -
+                     Math.Pow(Math.PI, 2) * Math.Pow(x, 2) * Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) +
+                     25 * Math.Pow(x, 2) + 10 * Math.PI * x * y * z - 50 * x * y * Math.Atan(y / x) +
+                     Math.Pow(Math.PI, 2) * Math.Pow(y, 4)) /
+              (Math.Pow(Math.PI, 2) * (Math.Pow(x, 4) + 2 * Math.Pow(x, 2) * Math.Pow(y, 2) + Math.Pow(y, 4)))
+            : 200 * (Math.Pow(Math.PI, 2) * Math.Pow(x, 4) +
+                     2 * Math.Pow(Math.PI, 2) * Math.Pow(x, 2) * Math.Pow(y, 2) -
+                     Math.Pow(Math.PI, 2) * Math.Pow(x, 2) * Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) +
+                     25 * Math.Pow(x, 2) + 10 * Math.PI * x * y * z - 50 * x * y * Math.Atan(y / x) -
+                     50 * Math.PI * x * y +
+                     Math.Pow(Math.PI, 2) * Math.Pow(y, 4)) / (Math.Pow(Math.PI, 2) *
+                                                               (Math.Pow(x, 4) + 2 * Math.Pow(x, 2) * Math.Pow(y, 2) +
+                                                                Math.Pow(y, 4)));
+        const double h22 = 202;
+        return [h00, h11, h22];
+    };
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/MinuitTutorialProblems/GoldsteinPriceProblem.cs b/test/ExampleProblems/MinuitTutorialProblems/GoldsteinPriceProblem.cs
new file mode 100644
index 0000000..cc898c3
--- /dev/null
+++ b/test/ExampleProblems/MinuitTutorialProblems/GoldsteinPriceProblem.cs
@@ -0,0 +1,117 @@
+namespace ExampleProblems.MinuitTutorialProblems;
+
+public class GoldsteinPriceProblem(DerivativeConfiguration derivativeConfiguration)
+    : MinuitTutorialProblem(
+        Parameters,
+        InitialValues,
+        OptimumValues,
+        Function,
+        Gradient,
+        Hessian,
+        HessianDiagonal,
+        derivativeConfiguration)
+{
+    // see section 7.5 in https://seal.web.cern.ch/seal/documents/minuit/mntutorial.pdf
+
+    private static readonly IReadOnlyList<string> Parameters = ["x", "y"];
+    private static readonly IReadOnlyList<double> InitialValues = [-0.4, -0.6];
+    private static readonly IReadOnlyList<double> OptimumValues = [0, -1];
+
+    private static readonly Func<IReadOnlyList<double>, double> Function = p =>
+    {
+        var (x, y) = (p[0], p[1]);
+        return (Math.Pow(2 * x - 3 * y, 2) *
+                   (12 * Math.Pow(x, 2) - 36 * x * y - 32 * x + 27 * Math.Pow(y, 2) + 48 * y + 18) + 30) *
+               (Math.Pow(x + y + 1, 2) * (3 * Math.Pow(x, 2) + 6 * x * y - 14 * x + 3 * Math.Pow(y, 2) - 14 * y + 19) +
+                1);
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> Gradient = p =>
+    {
+        var (x, y) = (p[0], p[1]);
+        var g0 =
+            4 * (2 * x - 3 * y) *
+            (Math.Pow(x + y + 1, 2) * (3 * Math.Pow(x, 2) + 6 * x * y - 14 * x + 3 * Math.Pow(y, 2) - 14 * y + 19) +
+             1) * (12 * Math.Pow(x, 2) - 36 * x * y - 32 * x + 27 * Math.Pow(y, 2) + 48 * y -
+                (2 * x - 3 * y) * (-6 * x + 9 * y + 8) + 18) + 2 *
+            (Math.Pow(2 * x - 3 * y, 2) *
+                (12 * Math.Pow(x, 2) - 36 * x * y - 32 * x + 27 * Math.Pow(y, 2) + 48 * y + 18) + 30) * (x + y + 1) *
+            (3 * Math.Pow(x, 2) + 6 * x * y - 14 * x + 3 * Math.Pow(y, 2) - 14 * y + (x + y + 1) * (3 * x + 3 * y - 7) +
+             19);
+        var g1 =
+            6 * (2 * x - 3 * y) *
+            (Math.Pow(x + y + 1, 2) * (3 * Math.Pow(x, 2) + 6 * x * y - 14 * x + 3 * Math.Pow(y, 2) - 14 * y + 19) +
+             1) * (-12 * Math.Pow(x, 2) + 36 * x * y + 32 * x - 27 * Math.Pow(y, 2) - 48 * y +
+                (2 * x - 3 * y) * (-6 * x + 9 * y + 8) - 18) + 2 *
+            (Math.Pow(2 * x - 3 * y, 2) *
+                (12 * Math.Pow(x, 2) - 36 * x * y - 32 * x + 27 * Math.Pow(y, 2) + 48 * y + 18) + 30) * (x + y + 1) *
+            (3 * Math.Pow(x, 2) + 6 * x * y - 14 * x + 3 * Math.Pow(y, 2) - 14 * y + (x + y + 1) * (3 * x + 3 * y - 7) +
+             19);
+        return [g0, g1];
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> Hessian = p =>
+    {
+        var (x, y) = (p[0], p[1]);
+        var h00 = 8064 * Math.Pow(x, 6) - 12096 * Math.Pow(x, 5) * y - 32256 * Math.Pow(x, 5) -
+                  19440 * Math.Pow(x, 4) * Math.Pow(y, 2) + 40320 * Math.Pow(x, 4) * y + 28560 * Math.Pow(x, 4) +
+                  24480 * Math.Pow(x, 3) * Math.Pow(y, 3) + 51840 * Math.Pow(x, 3) * Math.Pow(y, 2) -
+                  3360 * Math.Pow(x, 3) * y + 26880 * Math.Pow(x, 3) + 15660 * Math.Pow(x, 2) * Math.Pow(y, 4) -
+                  48960 * Math.Pow(x, 2) * Math.Pow(y, 3) - 64440 * Math.Pow(x, 2) * Math.Pow(y, 2) -
+                  92160 * Math.Pow(x, 2) * y - 29448 * Math.Pow(x, 2) - 11016 * x * Math.Pow(y, 5) -
+                  20880 * x * Math.Pow(y, 4) + 7440 * x * Math.Pow(y, 3) + 59040 * x * Math.Pow(y, 2) + 34704 * x * y -
+                  6432 * x - 2916 * Math.Pow(y, 6) + 7344 * Math.Pow(y, 5) + 17460 * Math.Pow(y, 4) +
+                  10080 * Math.Pow(y, 3) +
+                  15552 * Math.Pow(y, 2) + 14688 * y + 2520;
+        var h01 = -2016 * Math.Pow(x, 6) - 7776 * Math.Pow(x, 5) * y + 8064 * Math.Pow(x, 5) +
+                  18360 * Math.Pow(x, 4) * Math.Pow(y, 2) + 25920 * Math.Pow(x, 4) * y - 840 * Math.Pow(x, 4) +
+                  20880 * Math.Pow(x, 3) * Math.Pow(y, 3) - 48960 * Math.Pow(x, 3) * Math.Pow(y, 2) -
+                  42960 * Math.Pow(x, 3) * y - 30720 * Math.Pow(x, 3) - 27540 * Math.Pow(x, 2) * Math.Pow(y, 4) -
+                  41760 * Math.Pow(x, 2) * Math.Pow(y, 3) + 11160 * Math.Pow(x, 2) * Math.Pow(y, 2) +
+                  59040 * Math.Pow(x, 2) * y + 17352 * Math.Pow(x, 2) - 17496 * x * Math.Pow(y, 5) +
+                  36720 * x * Math.Pow(y, 4) + 69840 * x * Math.Pow(y, 3) + 30240 * x * Math.Pow(y, 2) + 31104 * x * y +
+                  14688 * x + 6804 * Math.Pow(y, 6) + 11664 * Math.Pow(y, 5) - 5940 * Math.Pow(y, 4) -
+                  47520 * Math.Pow(y, 3) - 70848 * Math.Pow(y, 2) - 38592 * y - 4680;
+        var h11 = -1296 * Math.Pow(x, 6) + 7344 * Math.Pow(x, 5) * y + 5184 * Math.Pow(x, 5) +
+                  15660 * Math.Pow(x, 4) * Math.Pow(y, 2) - 24480 * Math.Pow(x, 4) * y - 10740 * Math.Pow(x, 4) -
+                  36720 * Math.Pow(x, 3) * Math.Pow(y, 3) - 41760 * Math.Pow(x, 3) * Math.Pow(y, 2) +
+                  7440 * Math.Pow(x, 3) * y + 19680 * Math.Pow(x, 3) - 43740 * Math.Pow(x, 2) * Math.Pow(y, 4) +
+                  73440 * Math.Pow(x, 2) * Math.Pow(y, 3) + 104760 * Math.Pow(x, 2) * Math.Pow(y, 2) +
+                  30240 * Math.Pow(x, 2) * y + 15552 * Math.Pow(x, 2) + 40824 * x * Math.Pow(y, 5) +
+                  58320 * x * Math.Pow(y, 4) - 23760 * x * Math.Pow(y, 3) - 142560 * x * Math.Pow(y, 2) -
+                  141696 * x * y -
+                  38592 * x + 40824 * Math.Pow(y, 6) - 27216 * Math.Pow(y, 5) - 132840 * Math.Pow(y, 4) +
+                  38880 * Math.Pow(y, 3) + 172152 * Math.Pow(y, 2) + 73728 * y + 6120;
+        return
+        [
+            h00, h01,
+            h01, h11
+        ];
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> HessianDiagonal = p =>
+    {
+        var (x, y) = (p[0], p[1]);
+        var h00 = 8064 * Math.Pow(x, 6) - 12096 * Math.Pow(x, 5) * y - 32256 * Math.Pow(x, 5) -
+                  19440 * Math.Pow(x, 4) * Math.Pow(y, 2) + 40320 * Math.Pow(x, 4) * y + 28560 * Math.Pow(x, 4) +
+                  24480 * Math.Pow(x, 3) * Math.Pow(y, 3) + 51840 * Math.Pow(x, 3) * Math.Pow(y, 2) -
+                  3360 * Math.Pow(x, 3) * y + 26880 * Math.Pow(x, 3) + 15660 * Math.Pow(x, 2) * Math.Pow(y, 4) -
+                  48960 * Math.Pow(x, 2) * Math.Pow(y, 3) - 64440 * Math.Pow(x, 2) * Math.Pow(y, 2) -
+                  92160 * Math.Pow(x, 2) * y - 29448 * Math.Pow(x, 2) - 11016 * x * Math.Pow(y, 5) -
+                  20880 * x * Math.Pow(y, 4) + 7440 * x * Math.Pow(y, 3) + 59040 * x * Math.Pow(y, 2) + 34704 * x * y -
+                  6432 * x - 2916 * Math.Pow(y, 6) + 7344 * Math.Pow(y, 5) + 17460 * Math.Pow(y, 4) +
+                  10080 * Math.Pow(y, 3) +
+                  15552 * Math.Pow(y, 2) + 14688 * y + 2520;
+        var h11 = -1296 * Math.Pow(x, 6) + 7344 * Math.Pow(x, 5) * y + 5184 * Math.Pow(x, 5) +
+                  15660 * Math.Pow(x, 4) * Math.Pow(y, 2) - 24480 * Math.Pow(x, 4) * y - 10740 * Math.Pow(x, 4) -
+                  36720 * Math.Pow(x, 3) * Math.Pow(y, 3) - 41760 * Math.Pow(x, 3) * Math.Pow(y, 2) +
+                  7440 * Math.Pow(x, 3) * y + 19680 * Math.Pow(x, 3) - 43740 * Math.Pow(x, 2) * Math.Pow(y, 4) +
+                  73440 * Math.Pow(x, 2) * Math.Pow(y, 3) + 104760 * Math.Pow(x, 2) * Math.Pow(y, 2) +
+                  30240 * Math.Pow(x, 2) * y + 15552 * Math.Pow(x, 2) + 40824 * x * Math.Pow(y, 5) +
+                  58320 * x * Math.Pow(y, 4) - 23760 * x * Math.Pow(y, 3) - 142560 * x * Math.Pow(y, 2) -
+                  141696 * x * y -
+                  38592 * x + 40824 * Math.Pow(y, 6) - 27216 * Math.Pow(y, 5) - 132840 * Math.Pow(y, 4) +
+                  38880 * Math.Pow(y, 3) + 172152 * Math.Pow(y, 2) + 73728 * y + 6120;
+        return [h00, h11];
+    };
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/MinuitTutorialProblems/MinuitTutorialProblem.cs b/test/ExampleProblems/MinuitTutorialProblems/MinuitTutorialProblem.cs
new file mode 100644
index 0000000..f33aa4c
--- /dev/null
+++ b/test/ExampleProblems/MinuitTutorialProblems/MinuitTutorialProblem.cs
@@ -0,0 +1,32 @@
+using minuit2.net;
+using minuit2.net.CostFunctions;
+using static ExampleProblems.DerivativeConfiguration;
+using static minuit2.net.ParameterConfiguration;
+
+namespace ExampleProblems.MinuitTutorialProblems;
+
+public abstract class MinuitTutorialProblem(
+    IReadOnlyList<string> parameters,
+    IReadOnlyList<double> initialValues,
+    IReadOnlyList<double> optimumValues,
+    Func<IReadOnlyList<double>, double> function,
+    Func<IReadOnlyList<double>, IReadOnlyList<double>> gradient,
+    Func<IReadOnlyList<double>, IReadOnlyList<double>> hessian,
+    Func<IReadOnlyList<double>, IReadOnlyList<double>> hessianDiagonal,
+    DerivativeConfiguration derivativeConfiguration)
+    : IConfiguredProblem
+{
+    public ICostFunction Cost { get; } = derivativeConfiguration switch
+    {
+        WithoutDerivatives => new DirectCostFunction(parameters, function),
+        WithGradient => new DirectCostFunction(parameters, function, gradient),
+        WithGradientAndHessian => new DirectCostFunction(parameters, function, gradient, hessian),
+        WithGradientHessianAndHessianDiagonal => new DirectCostFunction(parameters, function, gradient, hessian, hessianDiagonal),
+        _ => throw new ArgumentOutOfRangeException(nameof(derivativeConfiguration), derivativeConfiguration, null)
+    };
+
+    public IReadOnlyCollection<double> OptimumParameterValues { get; } = optimumValues;
+
+    public IReadOnlyCollection<ParameterConfiguration> ParameterConfigurations { get; } =
+        parameters.Zip(initialValues, (name, value) => Variable(name, value)).ToList();
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/MinuitTutorialProblems/PowellProblem.cs b/test/ExampleProblems/MinuitTutorialProblems/PowellProblem.cs
new file mode 100644
index 0000000..b025d76
--- /dev/null
+++ b/test/ExampleProblems/MinuitTutorialProblems/PowellProblem.cs
@@ -0,0 +1,67 @@
+namespace ExampleProblems.MinuitTutorialProblems;
+
+public class PowellProblem(DerivativeConfiguration derivativeConfiguration)
+    : MinuitTutorialProblem(
+        Parameters,
+        InitialValues,
+        OptimumValues,
+        Function,
+        Gradient,
+        Hessian,
+        HessianDiagonal,
+        derivativeConfiguration)
+{
+    // see section 7.3 in https://seal.web.cern.ch/seal/documents/minuit/mntutorial.pdf
+
+    private static readonly IReadOnlyList<string> Parameters = ["w", "x", "y", "z"];
+    private static readonly IReadOnlyList<double> InitialValues = [3, -1, 0, 1];
+    private static readonly IReadOnlyList<double> OptimumValues = [0, 0, 0, 0];
+
+    private static readonly Func<IReadOnlyList<double>, double> Function = p =>
+    {
+        var (w, x, y, z) = (p[0], p[1], p[2], p[3]);
+        return Math.Pow(w + 10 * x, 2) + 10 * Math.Pow(w - z, 4) + Math.Pow(x - 2 * y, 4) + 5 * Math.Pow(y - z, 2);
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> Gradient = p =>
+    {
+        var (w, x, y, z) = (p[0], p[1], p[2], p[3]);
+        var g0 = 2 * w + 20 * x + 40 * Math.Pow(w - z, 3);
+        var g1 = 20 * w + 200 * x + 4 * Math.Pow(x - 2 * y, 3);
+        var g2 = 10 * y - 10 * z - 8 * Math.Pow(x - 2 * y, 3);
+        var g3 = -10 * y + 10 * z - 40 * Math.Pow(w - z, 3);
+        return [g0, g1, g2, g3];
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> Hessian = p =>
+    {
+        var (w, x, y, z) = (p[0], p[1], p[2], p[3]);
+        var h00 = 120 * Math.Pow(w - z, 2) + 2;
+        const double h01 = 20;
+        const double h02 = 0;
+        var h03 = -120 * Math.Pow(w - z, 2);
+        var h11 = 12 * Math.Pow(x - 2 * y, 2) + 200;
+        var h12 = -24 * Math.Pow(x - 2 * y, 2);
+        const double h13 = 0;
+        var h22 = 48 * Math.Pow(x - 2 * y, 2) + 10;
+        const double h23 = -10;
+        var h33 = 120 * Math.Pow(w - z, 2) + 10;
+        return
+        [
+            h00, h01, h02, h03,
+            h01, h11, h12, h13,
+            h02, h12, h22, h23,
+            h03, h13, h23, h33
+        ];
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> HessianDiagonal = p =>
+    {
+        var (w, x, y, z) = (p[0], p[1], p[2], p[3]);
+        var h00 = 120 * Math.Pow(w - z, 2) + 2;
+        var h11 = 12 * Math.Pow(x - 2 * y, 2) + 200;
+        var h22 = 48 * Math.Pow(x - 2 * y, 2) + 10;
+        var h33 = 120 * Math.Pow(w - z, 2) + 10;
+        return [h00, h11, h22, h33];
+    };
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/MinuitTutorialProblems/RosenbrockProblem.cs b/test/ExampleProblems/MinuitTutorialProblems/RosenbrockProblem.cs
new file mode 100644
index 0000000..f883a45
--- /dev/null
+++ b/test/ExampleProblems/MinuitTutorialProblems/RosenbrockProblem.cs
@@ -0,0 +1,54 @@
+namespace ExampleProblems.MinuitTutorialProblems;
+
+public class RosenbrockProblem(DerivativeConfiguration derivativeConfiguration)
+    : MinuitTutorialProblem(
+        Parameters,
+        InitialValues,
+        OptimumValues,
+        Function,
+        Gradient,
+        Hessian,
+        HessianDiagonal,
+        derivativeConfiguration)
+{
+    // see section 7.1 in https://seal.web.cern.ch/seal/documents/minuit/mntutorial.pdf
+
+    private static readonly IReadOnlyList<string> Parameters = ["x", "y"];
+    private static readonly IReadOnlyList<double> InitialValues = [-1.2, 1];
+    private static readonly IReadOnlyList<double> OptimumValues = [1, 1];
+
+    private static readonly Func<IReadOnlyList<double>, double> Function = p =>
+    {
+        var (x, y) = (p[0], p[1]);
+        return Math.Pow(1 - x, 2) + 100 * Math.Pow(-Math.Pow(x, 2) + y, 2);
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> Gradient = p =>
+    {
+        var (x, y) = (p[0], p[1]);
+        var g0 = 400 * x * (Math.Pow(x, 2) - y) + 2 * x - 2;
+        var g1 = -200 * Math.Pow(x, 2) + 200 * y;
+        return [g0, g1];
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> Hessian = p =>
+    {
+        var (x, y) = (p[0], p[1]);
+        var h00 = 1200 * Math.Pow(x, 2) - 400 * y + 2;
+        var h01 = -400 * x;
+        const double h11 = 200;
+        return
+        [
+            h00, h01,
+            h01, h11
+        ];
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> HessianDiagonal = p =>
+    {
+        var (x, y) = (p[0], p[1]);
+        var h00 = 1200 * Math.Pow(x, 2) - 400 * y + 2;
+        const double h11 = 200;
+        return [h00, h11];
+    };
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/MinuitTutorialProblems/WoodProblem.cs b/test/ExampleProblems/MinuitTutorialProblems/WoodProblem.cs
new file mode 100644
index 0000000..3ac49f9
--- /dev/null
+++ b/test/ExampleProblems/MinuitTutorialProblems/WoodProblem.cs
@@ -0,0 +1,69 @@
+namespace ExampleProblems.MinuitTutorialProblems;
+
+public class WoodProblem(DerivativeConfiguration derivativeConfiguration)
+    : MinuitTutorialProblem(
+        Parameters,
+        InitialValues,
+        OptimumValues,
+        Function,
+        Gradient,
+        Hessian,
+        HessianDiagonal,
+        derivativeConfiguration)
+{
+    // see section 7.2 in https://seal.web.cern.ch/seal/documents/minuit/mntutorial.pdf
+
+    private static readonly IReadOnlyList<string> Parameters = ["w", "x", "y", "z"];
+    private static readonly IReadOnlyList<double> InitialValues = [-3, -1, -3, -1];
+    private static readonly IReadOnlyList<double> OptimumValues = [1, 1, 1, 1];
+
+    private static readonly Func<IReadOnlyList<double>, double> Function = p =>
+    {
+        var (w, x, y, z) = (p[0], p[1], p[2], p[3]);
+        return Math.Pow(1 - y, 2) + Math.Pow(w - 1, 2) + 100 * Math.Pow(-Math.Pow(w, 2) + x, 2) +
+               10.1 * Math.Pow(x - 1, 2) + (19.8 * x - 19.8) * (z - 1) + 90 * Math.Pow(-Math.Pow(y, 2) + z, 2) +
+               10.1 * Math.Pow(z - 1, 2);
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> Gradient = p =>
+    {
+        var (w, x, y, z) = (p[0], p[1], p[2], p[3]);
+        var g0 = 400 * w * (Math.Pow(w, 2) - x) + 2 * w - 2;
+        var g1 = -200 * Math.Pow(w, 2) + 220.2 * x + 19.8 * z - 40.0;
+        var g2 = 360 * y * (Math.Pow(y, 2) - z) + 2 * y - 2;
+        var g3 = 19.8 * x - 180 * Math.Pow(y, 2) + 200.2 * z - 40.0;
+        return [g0, g1, g2, g3];
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> Hessian = p =>
+    {
+        var (w, x, y, z) = (p[0], p[1], p[2], p[3]);
+        var h00 = 1200 * Math.Pow(w, 2) - 400 * x + 2;
+        var h01 = -400 * w;
+        const double h02 = 0;
+        const double h03 = 0;
+        const double h11 = 220.2;
+        const double h12 = 0;
+        const double h13 = 19.8;
+        var h22 = 1080 * Math.Pow(y, 2) - 360 * z + 2;
+        var h23 = -360 * y;
+        const double h33 = 200.2;
+        return
+        [
+            h00, h01, h02, h03,
+            h01, h11, h12, h13,
+            h02, h12, h22, h23,
+            h03, h13, h23, h33
+        ];
+    };
+
+    private static readonly Func<IReadOnlyList<double>, IReadOnlyList<double>> HessianDiagonal = p =>
+    {
+        var (w, x, y, z) = (p[0], p[1], p[2], p[3]);
+        var h00 = 1200 * Math.Pow(w, 2) - 400 * x + 2;
+        const double h11 = 220.2;
+        var h22 = 1080 * Math.Pow(y, 2) - 360 * z + 2;
+        const double h33 = 200.2;
+        return [h00, h11, h22, h33];
+    };
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/NISTProblems/Bennett5Problem.cs b/test/ExampleProblems/NISTProblems/Bennett5Problem.cs
new file mode 100644
index 0000000..0285e8a
--- /dev/null
+++ b/test/ExampleProblems/NISTProblems/Bennett5Problem.cs
@@ -0,0 +1,104 @@
+namespace ExampleProblems.NISTProblems;
+
+public class Bennett5Problem(DerivativeConfiguration derivativeConfiguration)
+    : NistProblem(
+        X,
+        Y,
+        Parameters,
+        InitialValues,
+        OptimumValues,
+        Model,
+        ModelGradient,
+        ModelHessian,
+        ModelHessianDiagonal,
+        derivativeConfiguration)
+{
+    // see https://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml
+
+    private static readonly IReadOnlyList<double> X =
+    [
+        7.447168, 8.102586, 8.452547, 8.711278, 8.916774, 9.087155, 9.23259, 9.359535, 9.472166, 9.573384, 9.665293,
+        9.749461, 9.827092, 9.899128, 9.966321, 10.02928, 10.08851, 10.14443, 10.19738, 10.24767, 10.29556, 10.34125,
+        10.38495, 10.42682, 10.467, 10.50564, 10.54283, 10.57869, 10.61331, 10.64678, 10.67915, 10.71052, 10.74092,
+        10.77044, 10.7991, 10.82697, 10.85408, 10.88047, 10.90619, 10.93126, 10.95572, 10.97959, 11.00291, 11.0257,
+        11.04798, 11.06977, 11.0911, 11.11198, 11.13244, 11.15248, 11.17213, 11.19141, 11.21031, 11.22887, 11.24709,
+        11.26498, 11.28256, 11.29984, 11.31682, 11.33352, 11.34994, 11.3661, 11.382, 11.39766, 11.41307, 11.42824,
+        11.4432, 11.45793, 11.47244, 11.48675, 11.50086, 11.51477, 11.52849, 11.54202, 11.55538, 11.56855, 11.58156,
+        11.59442, 11.607121, 11.61964, 11.632, 11.64421, 11.65628, 11.6682, 11.67998, 11.69162, 11.70313, 11.71451,
+        11.72576, 11.73688, 11.74789, 11.75878, 11.76955, 11.7802, 11.79073, 11.80116, 11.81148, 11.8217, 11.83181,
+        11.84182, 11.85173, 11.86155, 11.87127, 11.88089, 11.89042, 11.89987, 11.90922, 11.91849, 11.92768, 11.93678,
+        11.94579, 11.95473, 11.96359, 11.97237, 11.98107, 11.9897, 11.99826, 12.00674, 12.01515, 12.02349, 12.03176,
+        12.03997, 12.0481, 12.05617, 12.06418, 12.07212, 12.08001, 12.08782, 12.09558, 12.10328, 12.11092, 12.1185,
+        12.12603, 12.1335, 12.14091, 12.14827, 12.15557, 12.16283, 12.17003, 12.17717, 12.18427, 12.19132, 12.19832,
+        12.20527, 12.21217, 12.21903, 12.22584, 12.2326, 12.23932, 12.24599, 12.25262, 12.2592, 12.26575, 12.27224
+    ];
+
+    private static readonly IReadOnlyList<double> Y =
+    [
+        -34.834702, -34.3932, -34.152901, -33.979099, -33.845901, -33.732899, -33.640301, -33.5592, -33.486801,
+        -33.4231, -33.365101, -33.313, -33.260899, -33.2174, -33.176899, -33.139198, -33.101601, -33.066799, -33.035,
+        -33.003101, -32.971298, -32.942299, -32.916302, -32.890202, -32.864101, -32.841, -32.817799, -32.797501,
+        -32.7743, -32.757, -32.733799, -32.7164, -32.6991, -32.678799, -32.6614, -32.644001, -32.626701, -32.612202,
+        -32.597698, -32.583199, -32.568699, -32.554298, -32.539799, -32.525299, -32.510799, -32.499199, -32.487598,
+        -32.473202, -32.461601, -32.435501, -32.435501, -32.4268, -32.4123, -32.400799, -32.392101, -32.380501,
+        -32.366001, -32.3573, -32.348598, -32.339901, -32.3284, -32.319698, -32.311001, -32.2994, -32.290699,
+        -32.282001, -32.2733, -32.264599, -32.256001, -32.247299, -32.238602, -32.2299, -32.224098, -32.215401,
+        -32.2038, -32.198002, -32.1894, -32.183601, -32.1749, -32.169102, -32.1633, -32.154598, -32.145901, -32.140099,
+        -32.131401, -32.125599, -32.119801, -32.111198, -32.1054, -32.096699, -32.0909, -32.088001, -32.0793,
+        -32.073502, -32.067699, -32.061901, -32.056099, -32.050301, -32.044498, -32.038799, -32.033001, -32.027199,
+        -32.0243, -32.018501, -32.012699, -32.004002, -32.001099, -31.9953, -31.9895, -31.9837, -31.9779, -31.972099,
+        -31.969299, -31.963501, -31.957701, -31.9519, -31.9461, -31.9403, -31.937401, -31.931601, -31.9258, -31.922899,
+        -31.917101, -31.911301, -31.9084, -31.902599, -31.8969, -31.893999, -31.888201, -31.8853, -31.882401, -31.8766,
+        -31.873699, -31.867901, -31.862101, -31.8592, -31.8563, -31.8505, -31.8447, -31.841801, -31.8389, -31.833099,
+        -31.8302, -31.827299, -31.8216, -31.818701, -31.812901, -31.809999, -31.8071, -31.8013, -31.798401, -31.7955,
+        -31.7897, -31.7868
+    ];
+
+    private static readonly IReadOnlyList<string> Parameters =
+        ["b1", "b2", "b3"];
+
+    private static readonly IReadOnlyList<double> InitialValues =
+        [-2000, 50, 0.8];
+
+    private static readonly IReadOnlyList<double> OptimumValues =
+        [-2.5235058043E+03, 4.6736564644E+01, 9.3218483193E-01];
+
+    private static readonly Func<double, IReadOnlyList<double>, double> Model = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        return b1 * Math.Pow(b2 + x, -1 / b3);
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelGradient = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        var g0 = Math.Pow(b2 + x, -1 / b3);
+        var g1 = -b1 * Math.Pow(b2 + x, -(b3 + 1) / b3) / b3;
+        var g2 = b1 * Math.Pow(b2 + x, -1 / b3) * Math.Log(b2 + x) / Math.Pow(b3, 2);
+        return [g0, g1, g2];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessian = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        var h01 = -Math.Pow(b2 + x, -(b3 + 1) / b3) / b3;
+        var h02 = Math.Pow(b2 + x, -1 / b3) * Math.Log(b2 + x) / Math.Pow(b3, 2);
+        var h11 = Math.Pow(b2 + x, -2 - 1 / b3) * (b1 * b3 + b1) / Math.Pow(b3, 2);
+        var h12 = b1 * Math.Pow(b2 + x, -(b3 + 1) / b3) * (b3 - Math.Log(b2 + x)) / Math.Pow(b3, 3);
+        var h22 = b1 * Math.Pow(b2 + x, -1 / b3) * (-2 * b3 + Math.Log(b2 + x)) * Math.Log(b2 + x) / Math.Pow(b3, 4);
+        return
+        [
+            0.0, h01, h02,
+            h01, h11, h12,
+            h02, h12, h22
+        ];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessianDiagonal = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        var h11 = Math.Pow(b2 + x, -2 - 1 / b3) * (b1 * b3 + b1) / Math.Pow(b3, 2);
+        var h22 = b1 * Math.Pow(b2 + x, -1 / b3) * (-2 * b3 + Math.Log(b2 + x)) * Math.Log(b2 + x) / Math.Pow(b3, 4);
+        return [0, h11, h22];
+    };
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/NISTProblems/BoxBodProblem.cs b/test/ExampleProblems/NISTProblems/BoxBodProblem.cs
new file mode 100644
index 0000000..48b039e
--- /dev/null
+++ b/test/ExampleProblems/NISTProblems/BoxBodProblem.cs
@@ -0,0 +1,56 @@
+namespace ExampleProblems.NISTProblems;
+
+public class BoxBodProblem(DerivativeConfiguration derivativeConfiguration)
+    : NistProblem(
+        X,
+        Y,
+        Parameters,
+        InitialValues,
+        OptimumValues,
+        Model,
+        ModelGradient,
+        ModelHessian,
+        ModelHessianDiagonal,
+        derivativeConfiguration)
+{
+    // see https://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml
+
+    private static readonly IReadOnlyList<double> X = [1, 2, 3, 5, 7, 10];
+    private static readonly IReadOnlyList<double> Y = [109, 149, 149, 191, 213, 224];
+    private static readonly IReadOnlyList<string> Parameters = ["b1", "b2"];
+    private static readonly IReadOnlyList<double> InitialValues = [1, 1];
+    private static readonly IReadOnlyList<double> OptimumValues = [2.1380940889E+02, 5.4723748542E-01];
+
+    private static readonly Func<double, IReadOnlyList<double>, double> Model = (x, p) =>
+    {
+        var (b1, b2) = (p[0], p[1]);
+        return b1 * (1 - Math.Exp(-b2 * x));
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelGradient = (x, p) =>
+    {
+        var (b1, b2) = (p[0], p[1]);
+        var g0 = 1 - Math.Exp(-b2 * x);
+        var g1 = b1 * x * Math.Exp(-b2 * x);
+        return [g0, g1];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessian = (x, p) =>
+    {
+        var (b1, b2) = (p[0], p[1]);
+        var h01 = x * Math.Exp(-b2 * x);
+        var h11 = -b1 * Math.Pow(x, 2) * Math.Exp(-b2 * x);
+        return
+        [
+            0.0, h01,
+            h01, h11
+        ];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessianDiagonal = (x, p) =>
+    {
+        var (b1, b2) = (p[0], p[1]);
+        var h11 = -b1 * Math.Pow(x, 2) * Math.Exp(-b2 * x);
+        return [0, h11];
+    };
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/NISTProblems/Eckerle4Problem.cs b/test/ExampleProblems/NISTProblems/Eckerle4Problem.cs
new file mode 100644
index 0000000..f4a92f2
--- /dev/null
+++ b/test/ExampleProblems/NISTProblems/Eckerle4Problem.cs
@@ -0,0 +1,90 @@
+using System.Diagnostics.CodeAnalysis;
+
+namespace ExampleProblems.NISTProblems;
+
+[SuppressMessage("ReSharper", "UnusedType.Global", Justification = "Crashes the native minimization processes. Needs further investigation.")]
+public class Eckerle4Problem(DerivativeConfiguration derivativeConfiguration)
+    : NistProblem(
+        X,
+        Y,
+        Parameters,
+        InitialValues,
+        OptimumValues,
+        Model,
+        ModelGradient,
+        ModelHessian,
+        ModelHessianDiagonal,
+        derivativeConfiguration)
+{
+    // see https://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml
+
+    private static readonly IReadOnlyList<double> X =
+    [
+        400.0, 405.0, 410.0, 415.0, 420.0, 425.0, 430.0, 435.0, 436.5, 438.0, 439.5, 441.0, 442.5, 444.0, 445.5, 447.0,
+        448.5, 450.0, 451.5, 453.0, 454.5, 456.0, 457.5, 459.0, 460.5, 462.0, 463.5, 465.0, 470.0, 475.0, 480.0, 485.0,
+        490.0, 495.0, 500.0
+    ];
+
+    private static readonly IReadOnlyList<double> Y =
+    [
+        0.0001575, 0.0001699, 0.000235, 0.0003102, 0.0004917, 0.000871, 0.0017418, 0.00464, 0.0065895, 0.0097302,
+        0.0149002, 0.023731, 0.0401683, 0.0712559, 0.1264458, 0.2073413, 0.2902366, 0.3445623, 0.3698049, 0.3668534,
+        0.3106727, 0.2078154, 0.1164354, 0.0616764, 0.03372, 0.0194023, 0.0117831, 0.0074357, 0.0022732, 0.00088,
+        0.0004579, 0.0002345, 0.0001586, 0.0001143, 7.1e-05
+    ];
+
+    private static readonly IReadOnlyList<string> Parameters =
+        ["b1", "b2", "b3"];
+
+    private static readonly IReadOnlyList<double> InitialValues =
+        [1, 10, 500];
+
+    private static readonly IReadOnlyList<double> OptimumValues =
+        [1.5543827178E+00, 4.0888321754E+00, 4.5154121844E+02];
+
+    private static readonly Func<double, IReadOnlyList<double>, double> Model = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        return b1 * Math.Exp(-0.5 * Math.Pow(-b3 + x, 2) / Math.Pow(b2, 2)) / b2;
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelGradient = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        var g0 = Math.Exp(-0.5 * Math.Pow(b3 - x, 2) / Math.Pow(b2, 2)) / b2;
+        var g1 = b1 * (-Math.Pow(b2, 2) + 1.0 * Math.Pow(b3 - x, 2)) *
+            Math.Exp(-0.5 * Math.Pow(b3 - x, 2) / Math.Pow(b2, 2)) / Math.Pow(b2, 4);
+        var g2 = -1.0 * b1 * (b3 - x) * Math.Exp(-0.5 * Math.Pow(b3 - x, 2) / Math.Pow(b2, 2)) / Math.Pow(b2, 3);
+        return [g0, g1, g2];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessian = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        var h01 = (-Math.Pow(b2, 2) + 1.0 * Math.Pow(b3 - x, 2)) *
+            Math.Exp(-0.5 * Math.Pow(b3 - x, 2) / Math.Pow(b2, 2)) / Math.Pow(b2, 4);
+        var h02 = 1.0 * (-b3 + x) * Math.Exp(-0.5 * Math.Pow(b3 - x, 2) / Math.Pow(b2, 2)) / Math.Pow(b2, 3);
+        var h11 = b1 * (2 * Math.Pow(b2, 4) - 5.0 * Math.Pow(b2, 2) * Math.Pow(b3 - x, 2) + 1.0 * Math.Pow(b3 - x, 4)) *
+            Math.Exp(-0.5 * Math.Pow(b3 - x, 2) / Math.Pow(b2, 2)) / Math.Pow(b2, 7);
+        var h12 = b1 * (3.0 * Math.Pow(b2, 2) - 1.0 * Math.Pow(b3 - x, 2)) * (b3 - x) *
+            Math.Exp(-0.5 * Math.Pow(b3 - x, 2) / Math.Pow(b2, 2)) / Math.Pow(b2, 6);
+        var h22 = 1.0 * b1 * (-Math.Pow(b2, 2) + Math.Pow(b3 - x, 2)) *
+            Math.Exp(-0.5 * Math.Pow(b3 - x, 2) / Math.Pow(b2, 2)) / Math.Pow(b2, 5);
+        return
+        [
+            0.0, h01, h02,
+            h01, h11, h12,
+            h02, h12, h22
+        ];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessianDiagonal = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        var h11 = b1 * (2 * Math.Pow(b2, 4) - 5.0 * Math.Pow(b2, 2) * Math.Pow(b3 - x, 2) + 1.0 * Math.Pow(b3 - x, 4)) *
+            Math.Exp(-0.5 * Math.Pow(b3 - x, 2) / Math.Pow(b2, 2)) / Math.Pow(b2, 7);
+        var h22 = 1.0 * b1 * (-Math.Pow(b2, 2) + Math.Pow(b3 - x, 2)) *
+            Math.Exp(-0.5 * Math.Pow(b3 - x, 2) / Math.Pow(b2, 2)) / Math.Pow(b2, 5);
+        return [0, h11, h22];
+    };
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/NISTProblems/Mgh09Problem.cs b/test/ExampleProblems/NISTProblems/Mgh09Problem.cs
new file mode 100644
index 0000000..ed5fb50
--- /dev/null
+++ b/test/ExampleProblems/NISTProblems/Mgh09Problem.cs
@@ -0,0 +1,76 @@
+namespace ExampleProblems.NISTProblems;
+
+public class Mgh09Problem(DerivativeConfiguration derivativeConfiguration)
+    : NistProblem(
+        X,
+        Y,
+        Parameters,
+        InitialValues,
+        OptimumValues,
+        Model,
+        ModelGradient,
+        ModelHessian,
+        ModelHessianDiagonal,
+        derivativeConfiguration)
+{
+    // see https://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml
+
+    private static readonly IReadOnlyList<double> X =
+        [4.0, 2.0, 1.0, 0.5, 0.25, 0.167, 0.125, 0.1, 0.0833, 0.0714, 0.0625];
+
+    private static readonly IReadOnlyList<double> Y =
+        [0.1957, 0.1947, 0.1735, 0.16, 0.0844, 0.0627, 0.0456, 0.0342, 0.0323, 0.0235, 0.0246];
+
+    private static readonly IReadOnlyList<string> Parameters =
+        ["b1", "b2", "b3", "b4"];
+
+    private static readonly IReadOnlyList<double> InitialValues =
+        [25, 39, 41.5, 39];
+
+    private static readonly IReadOnlyList<double> OptimumValues =
+        [1.9280693458E-01, 1.9128232873E-01, 1.2305650693E-01, 1.3606233068E-01];
+
+    private static readonly Func<double, IReadOnlyList<double>, double> Model = (x, p) =>
+    {
+        var (b1, b2, b3, b4) = (p[0], p[1], p[2], p[3]);
+        return b1 * (b2 * x + Math.Pow(x, 2)) / (b3 * x + b4 + Math.Pow(x, 2));
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelGradient = (x, p) =>
+    {
+        var (b1, b2, b3, b4) = (p[0], p[1], p[2], p[3]);
+        var g0 = x * (b2 + x) / (b3 * x + b4 + Math.Pow(x, 2));
+        var g1 = b1 * x / (b3 * x + b4 + Math.Pow(x, 2));
+        var g2 = -b1 * Math.Pow(x, 2) * (b2 + x) / Math.Pow(b3 * x + b4 + Math.Pow(x, 2), 2);
+        var g3 = -b1 * x * (b2 + x) / Math.Pow(b3 * x + b4 + Math.Pow(x, 2), 2);
+        return [g0, g1, g2, g3];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessian = (x, p) =>
+    {
+        var (b1, b2, b3, b4) = (p[0], p[1], p[2], p[3]);
+        var h01 = x / (b3 * x + b4 + Math.Pow(x, 2));
+        var h02 = -Math.Pow(x, 2) * (b2 + x) / Math.Pow(b3 * x + b4 + Math.Pow(x, 2), 2);
+        var h03 = x * (-b2 - x) / Math.Pow(b3 * x + b4 + Math.Pow(x, 2), 2);
+        var h12 = -b1 * Math.Pow(x, 2) / Math.Pow(b3 * x + b4 + Math.Pow(x, 2), 2);
+        var h13 = -b1 * x / Math.Pow(b3 * x + b4 + Math.Pow(x, 2), 2);
+        var h22 = 2 * b1 * Math.Pow(x, 3) * (b2 + x) / Math.Pow(b3 * x + b4 + Math.Pow(x, 2), 3);
+        var h23 = 2 * b1 * Math.Pow(x, 2) * (b2 + x) / Math.Pow(b3 * x + b4 + Math.Pow(x, 2), 3);
+        var h33 = 2 * b1 * x * (b2 + x) / Math.Pow(b3 * x + b4 + Math.Pow(x, 2), 3);
+        return
+        [
+            0.0, h01, h02, h03,
+            h01, 0.0, h12, h13,
+            h02, h12, h22, h23,
+            h03, h13, h23, h33
+        ];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessianDiagonal = (x, p) =>
+    {
+        var (b1, b2, b3, b4) = (p[0], p[1], p[2], p[3]);
+        var h22 = 2 * b1 * Math.Pow(x, 3) * (b2 + x) / Math.Pow(b3 * x + b4 + Math.Pow(x, 2), 3);
+        var h33 = 2 * b1 * x * (b2 + x) / Math.Pow(b3 * x + b4 + Math.Pow(x, 2), 3);
+        return [0, 0, h22, h33];
+    };
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/NISTProblems/Mgh10Problem.cs b/test/ExampleProblems/NISTProblems/Mgh10Problem.cs
new file mode 100644
index 0000000..0044d73
--- /dev/null
+++ b/test/ExampleProblems/NISTProblems/Mgh10Problem.cs
@@ -0,0 +1,77 @@
+using System.Diagnostics.CodeAnalysis;
+
+namespace ExampleProblems.NISTProblems;
+
+[SuppressMessage("ReSharper", "UnusedType.Global", Justification = "Crashes the native minimization processes. Needs further investigation.")]
+public class Mgh10Problem(DerivativeConfiguration derivativeConfiguration)
+    : NistProblem(
+        X,
+        Y,
+        Parameters,
+        InitialValues,
+        OptimumValues,
+        Model,
+        ModelGradient,
+        ModelHessian,
+        ModelHessianDiagonal,
+        derivativeConfiguration)
+{
+    // see https://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml
+
+    private static readonly IReadOnlyList<double> X =
+        [50.0, 55.0, 60.0, 65.0, 70.0, 75.0, 80.0, 85.0, 90.0, 95.0, 100.0, 105.0, 110.0, 115.0, 120.0, 125.0];
+
+    private static readonly IReadOnlyList<double> Y =
+    [
+        34780.0, 28610.0, 23650.0, 19630.0, 16370.0, 13720.0, 11540.0, 9744.0, 8261.0, 7030.0, 6005.0, 5147.0, 4427.0,
+        3820.0, 3307.0, 2872.0
+    ];
+
+    private static readonly IReadOnlyList<string> Parameters = 
+        ["b1", "b2", "b3"];
+    
+    private static readonly IReadOnlyList<double> InitialValues = 
+        [2, 400000, 25000];
+    
+    private static readonly IReadOnlyList<double> OptimumValues = 
+        [5.6096364710E-03, 6.1813463463E+03, 3.4522363462E+02];
+
+    private static readonly Func<double, IReadOnlyList<double>, double> Model = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        return b1 * Math.Exp(b2 / (b3 + x));
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelGradient = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        var g0 = Math.Exp(b2 / (b3 + x));
+        var g1 = b1 * Math.Exp(b2 / (b3 + x)) / (b3 + x);
+        var g2 = -b1 * b2 * Math.Exp(b2 / (b3 + x)) / Math.Pow(b3 + x, 2);
+        return [g0, g1, g2];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessian = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        var h01 = Math.Exp(b2 / (b3 + x)) / (b3 + x);
+        var h02 = -b2 * Math.Exp(b2 / (b3 + x)) / Math.Pow(b3 + x, 2);
+        var h11 = b1 * Math.Exp(b2 / (b3 + x)) / Math.Pow(b3 + x, 2);
+        var h12 = b1 * (-b2 - b3 - x) * Math.Exp(b2 / (b3 + x)) / Math.Pow(b3 + x, 3);
+        var h22 = b1 * b2 * (b2 + 2 * b3 + 2 * x) * Math.Exp(b2 / (b3 + x)) / Math.Pow(b3 + x, 4);
+        return
+        [
+            0.0, h01, h02,
+            h01, h11, h12,
+            h02, h12, h22
+        ];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessianDiagonal = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        var h11 = b1 * Math.Exp(b2 / (b3 + x)) / Math.Pow(b3 + x, 2);
+        var h22 = b1 * b2 * (b2 + 2 * b3 + 2 * x) * Math.Exp(b2 / (b3 + x)) / Math.Pow(b3 + x, 4);
+        return [0, h11, h22];
+    };
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/NISTProblems/NistProblem.cs b/test/ExampleProblems/NISTProblems/NistProblem.cs
new file mode 100644
index 0000000..1a511f1
--- /dev/null
+++ b/test/ExampleProblems/NISTProblems/NistProblem.cs
@@ -0,0 +1,34 @@
+using minuit2.net;
+using minuit2.net.CostFunctions;
+using static ExampleProblems.DerivativeConfiguration;
+using static minuit2.net.CostFunctions.CostFunction;
+
+namespace ExampleProblems.NISTProblems;
+
+public abstract class NistProblem(
+    IReadOnlyList<double> x,
+    IReadOnlyList<double> y,
+    IReadOnlyList<string> parameters,
+    IReadOnlyList<double> initialValues,
+    IReadOnlyList<double> optimumValues,
+    Func<double, IReadOnlyList<double>, double> model,
+    Func<double, IReadOnlyList<double>, IReadOnlyList<double>> modelGradient,
+    Func<double, IReadOnlyList<double>, IReadOnlyList<double>> modelHessian,
+    Func<double, IReadOnlyList<double>, IReadOnlyList<double>> modelHessianDiagonal,
+    DerivativeConfiguration modelDerivativeConfiguration)
+    : IConfiguredProblem
+{
+    public ICostFunction Cost { get; } = modelDerivativeConfiguration switch
+    {
+        WithoutDerivatives => LeastSquares(x, y, parameters, model),
+        WithGradient => LeastSquares(x, y, parameters, model, modelGradient),
+        WithGradientAndHessian => LeastSquares(x, y, parameters, model, modelGradient, modelHessian),
+        WithGradientHessianAndHessianDiagonal => LeastSquares(x, y, parameters, model, modelGradient, modelHessian, modelHessianDiagonal),
+        _ => throw new ArgumentOutOfRangeException(nameof(modelDerivativeConfiguration), modelDerivativeConfiguration, null)
+    };
+
+    public IReadOnlyCollection<double> OptimumParameterValues { get; } = optimumValues;
+
+    public IReadOnlyCollection<ParameterConfiguration> ParameterConfigurations { get; } =
+        parameters.Zip(initialValues, (name, value) => ParameterConfiguration.Variable(name, value)).ToList();
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/NISTProblems/Rat42Problem.cs b/test/ExampleProblems/NISTProblems/Rat42Problem.cs
new file mode 100644
index 0000000..f902fb5
--- /dev/null
+++ b/test/ExampleProblems/NISTProblems/Rat42Problem.cs
@@ -0,0 +1,69 @@
+namespace ExampleProblems.NISTProblems;
+
+public class Rat42Problem(DerivativeConfiguration derivativeConfiguration)
+    : NistProblem(
+        X,
+        Y,
+        Parameters,
+        InitialValues,
+        OptimumValues,
+        Model,
+        ModelGradient,
+        ModelHessian,
+        ModelHessianDiagonal,
+        derivativeConfiguration)
+{
+    // see https://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml
+
+    private static readonly IReadOnlyList<double> X = [9.0, 14.0, 21.0, 28.0, 42.0, 57.0, 63.0, 70.0, 79.0];
+    private static readonly IReadOnlyList<double> Y = [8.93, 10.8, 18.59, 22.33, 39.35, 56.11, 61.73, 64.62, 67.08];
+    private static readonly IReadOnlyList<string> Parameters = ["b1", "b2", "b3"];
+    private static readonly IReadOnlyList<double> InitialValues = [100, 1, 0.1];
+    private static readonly IReadOnlyList<double> OptimumValues = [7.2462237576E+01, 2.6180768402E+00, 6.7359200066E-02];
+
+    private static readonly Func<double, IReadOnlyList<double>, double> Model = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        return b1 / (Math.Exp(b2 - b3 * x) + 1);
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelGradient = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        var g0 = 1d / (Math.Exp(b2 - b3 * x) + 1);
+        var g1 = -1d / 4 * b1 / Math.Pow(Math.Cosh(1d / 2 * b2 - 1d / 2 * b3 * x), 2);
+        var g2 = 1d / 4 * b1 * x / Math.Pow(Math.Cosh(1d / 2 * b2 - 1d / 2 * b3 * x), 2);
+        return [g0, g1, g2];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessian = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        var h01 = -(1d / 4) / Math.Pow(Math.Cosh(1d / 2 * b2 - 1d / 2 * b3 * x), 2);
+        var h02 = 1d / 4 * x / Math.Pow(Math.Cosh(1d / 2 * b2 - 1d / 2 * b3 * x), 2);
+        var h11 = b1 * (-(Math.Exp(b2 - b3 * x) + 1) * Math.Exp(b2 - b3 * x) + 2 * Math.Exp(2 * b2 - 2 * b3 * x)) /
+                  Math.Pow(Math.Exp(b2 - b3 * x) + 1, 3);
+        var h12 = b1 * x * ((Math.Exp(b2 - b3 * x) + 1) * Math.Exp(b2 - b3 * x) - 2 * Math.Exp(2 * b2 - 2 * b3 * x)) /
+                  Math.Pow(Math.Exp(b2 - b3 * x) + 1, 3);
+        var h22 = b1 * Math.Pow(x, 2) *
+                  (-(Math.Exp(b2 - b3 * x) + 1) * Math.Exp(b2 - b3 * x) + 2 * Math.Exp(2 * b2 - 2 * b3 * x)) /
+                  Math.Pow(Math.Exp(b2 - b3 * x) + 1, 3);
+        return
+        [
+            0.0, h01, h02,
+            h01, h11, h12,
+            h02, h12, h22
+        ];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessianDiagonal = (x, p) =>
+    {
+        var (b1, b2, b3) = (p[0], p[1], p[2]);
+        var h11 = b1 * (-(Math.Exp(b2 - b3 * x) + 1) * Math.Exp(b2 - b3 * x) + 2 * Math.Exp(2 * b2 - 2 * b3 * x)) /
+                  Math.Pow(Math.Exp(b2 - b3 * x) + 1, 3);
+        var h22 = b1 * Math.Pow(x, 2) *
+                  (-(Math.Exp(b2 - b3 * x) + 1) * Math.Exp(b2 - b3 * x) + 2 * Math.Exp(2 * b2 - 2 * b3 * x)) /
+                  Math.Pow(Math.Exp(b2 - b3 * x) + 1, 3);
+        return [0, h11, h22];
+    };
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/NISTProblems/Rat43Problem.cs b/test/ExampleProblems/NISTProblems/Rat43Problem.cs
new file mode 100644
index 0000000..1342c55
--- /dev/null
+++ b/test/ExampleProblems/NISTProblems/Rat43Problem.cs
@@ -0,0 +1,100 @@
+namespace ExampleProblems.NISTProblems;
+
+public class Rat43Problem(DerivativeConfiguration derivativeConfiguration)
+    : NistProblem(
+        X,
+        Y,
+        Parameters,
+        InitialValues,
+        OptimumValues,
+        Model,
+        ModelGradient,
+        ModelHessian,
+        ModelHessianDiagonal,
+        derivativeConfiguration)
+{
+    // see https://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml
+
+    private static readonly IReadOnlyList<double> X =
+        [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0];
+
+    private static readonly IReadOnlyList<double> Y = 
+        [16.08, 33.83, 65.8, 97.2, 191.55, 326.2, 386.87, 520.53, 590.03, 651.92, 724.93, 699.56, 689.96, 637.56, 717.41];
+
+    private static readonly IReadOnlyList<string> Parameters =
+        ["b1", "b2", "b3", "b4"];
+
+    private static readonly IReadOnlyList<double> InitialValues =
+        [100, 10, 1, 1];
+
+    private static readonly IReadOnlyList<double> OptimumValues =
+        [6.9964151270E+02, 5.2771253025E+00, 7.5962938329E-01, 1.2792483859E+00];
+
+    private static readonly Func<double, IReadOnlyList<double>, double> Model = (x, p) =>
+    {
+        var (b1, b2, b3, b4) = (p[0], p[1], p[2], p[3]);
+        return b1 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, -1 / b4);
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelGradient = (x, p) =>
+    {
+        var (b1, b2, b3, b4) = (p[0], p[1], p[2], p[3]);
+        var g0 = Math.Pow(Math.Exp(b2 - b3 * x) + 1, -1 / b4);
+        var g1 = -b1 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, -(b4 + 1) / b4) * Math.Exp(b2 - b3 * x) / b4;
+        var g2 = b1 * x * Math.Pow(Math.Exp(b2 - b3 * x) + 1, -(b4 + 1) / b4) * Math.Exp(b2 - b3 * x) / b4;
+        var g3 = b1 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, -1 / b4) * Math.Log(Math.Exp(b2 - b3 * x) + 1) /
+                 Math.Pow(b4, 2);
+        return [g0, g1, g2, g3];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessian = (x, p) =>
+    {
+        var (b1, b2, b3, b4) = (p[0], p[1], p[2], p[3]);
+        var h01 = -Math.Pow(Math.Exp(b2 - b3 * x) + 1, -(b4 + 1) / b4) * Math.Exp(b2 - b3 * x) / b4;
+        var h02 = x * Math.Pow(Math.Exp(b2 - b3 * x) + 1, -(b4 + 1) / b4) * Math.Exp(b2 - b3 * x) / b4;
+        var h03 = Math.Pow(Math.Exp(b2 - b3 * x) + 1, -1 / b4) * Math.Log(Math.Exp(b2 - b3 * x) + 1) / Math.Pow(b4, 2);
+        var h11 = b1 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, -(5 * b4 + 3) / b4) *
+            (-b4 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, 2 * (2 * b4 + 1) / b4) * Math.Exp(b2 - b3 * x) +
+             b4 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, (3 * b4 + 2) / b4) * Math.Exp(2 * b2 - 2 * b3 * x) +
+             Math.Pow(Math.Exp(b2 - b3 * x) + 1, (3 * b4 + 2) / b4) * Math.Exp(2 * b2 - 2 * b3 * x)) / Math.Pow(b4, 2);
+        var h12 = b1 * x * Math.Pow(Math.Exp(b2 - b3 * x) + 1, -(5 * b4 + 3) / b4) *
+            (b4 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, 2 * (2 * b4 + 1) / b4) * Math.Exp(b2 - b3 * x) -
+             b4 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, (3 * b4 + 2) / b4) * Math.Exp(2 * b2 - 2 * b3 * x) -
+             Math.Pow(Math.Exp(b2 - b3 * x) + 1, (3 * b4 + 2) / b4) * Math.Exp(2 * b2 - 2 * b3 * x)) / Math.Pow(b4, 2);
+        var h13 = b1 * Math.Pow((Math.Exp(b2) + Math.Exp(b3 * x)) * Math.Exp(-b3 * x), -(b4 + 1) / b4) *
+                  (b4 - Math.Log((Math.Exp(b2) + Math.Exp(b3 * x)) * Math.Exp(-b3 * x))) * Math.Exp(b2 - b3 * x) /
+                  Math.Pow(b4, 3);
+        var h22 = b1 * Math.Pow(x, 2) * Math.Pow(Math.Exp(b2 - b3 * x) + 1, -(5 * b4 + 3) / b4) *
+            (-b4 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, 2 * (2 * b4 + 1) / b4) * Math.Exp(b2 - b3 * x) +
+             b4 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, (3 * b4 + 2) / b4) * Math.Exp(2 * b2 - 2 * b3 * x) +
+             Math.Pow(Math.Exp(b2 - b3 * x) + 1, (3 * b4 + 2) / b4) * Math.Exp(2 * b2 - 2 * b3 * x)) / Math.Pow(b4, 2);
+        var h23 = b1 * x * Math.Pow((Math.Exp(b2) + Math.Exp(b3 * x)) * Math.Exp(-b3 * x), -(b4 + 1) / b4) *
+                  (-b4 + Math.Log((Math.Exp(b2) + Math.Exp(b3 * x)) * Math.Exp(-b3 * x))) * Math.Exp(b2 - b3 * x) /
+                  Math.Pow(b4, 3);
+        var h33 = b1 * (-2 * b4 + Math.Log(Math.Exp(b2 - b3 * x) + 1)) * Math.Pow(Math.Exp(b2 - b3 * x) + 1, -1 / b4) *
+            Math.Log(Math.Exp(b2 - b3 * x) + 1) / Math.Pow(b4, 4);
+        return
+        [
+            0.0, h01, h02, h03,
+            h01, h11, h12, h13,
+            h02, h12, h22, h23,
+            h03, h13, h23, h33
+        ];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessianDiagonal = (x, p) =>
+    {
+        var (b1, b2, b3, b4) = (p[0], p[1], p[2], p[3]);
+        var h11 = b1 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, -(5 * b4 + 3) / b4) *
+            (-b4 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, 2 * (2 * b4 + 1) / b4) * Math.Exp(b2 - b3 * x) +
+             b4 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, (3 * b4 + 2) / b4) * Math.Exp(2 * b2 - 2 * b3 * x) +
+             Math.Pow(Math.Exp(b2 - b3 * x) + 1, (3 * b4 + 2) / b4) * Math.Exp(2 * b2 - 2 * b3 * x)) / Math.Pow(b4, 2);
+        var h22 = b1 * Math.Pow(x, 2) * Math.Pow(Math.Exp(b2 - b3 * x) + 1, -(5 * b4 + 3) / b4) *
+            (-b4 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, 2 * (2 * b4 + 1) / b4) * Math.Exp(b2 - b3 * x) +
+             b4 * Math.Pow(Math.Exp(b2 - b3 * x) + 1, (3 * b4 + 2) / b4) * Math.Exp(2 * b2 - 2 * b3 * x) +
+             Math.Pow(Math.Exp(b2 - b3 * x) + 1, (3 * b4 + 2) / b4) * Math.Exp(2 * b2 - 2 * b3 * x)) / Math.Pow(b4, 2);
+        var h33 = b1 * (-2 * b4 + Math.Log(Math.Exp(b2 - b3 * x) + 1)) * Math.Pow(Math.Exp(b2 - b3 * x) + 1, -1 / b4) *
+            Math.Log(Math.Exp(b2 - b3 * x) + 1) / Math.Pow(b4, 4);
+        return [0, h11, h22, h33];
+    };
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/NISTProblems/ThurberProblem.cs b/test/ExampleProblems/NISTProblems/ThurberProblem.cs
new file mode 100644
index 0000000..1115d7e
--- /dev/null
+++ b/test/ExampleProblems/NISTProblems/ThurberProblem.cs
@@ -0,0 +1,118 @@
+using System.Diagnostics.CodeAnalysis;
+
+namespace ExampleProblems.NISTProblems;
+
+public class ThurberProblem(DerivativeConfiguration derivativeConfiguration)
+    : NistProblem(
+        X,
+        Y,
+        Parameters,
+        InitialValues,
+        OptimumValues,
+        Model,
+        ModelGradient,
+        ModelHessian,
+        ModelHessianDiagonal,
+        derivativeConfiguration)
+{
+    // see https://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml
+
+    private static readonly IReadOnlyList<double> X =
+    [
+        -3.067, -2.981, -2.921, -2.912, -2.84, -2.797, -2.702, -2.699, -2.633, -2.481, -2.363, -2.322, -1.501, -1.46,
+        -1.274, -1.212, -1.1, -1.046, -0.915, -0.714, -0.566, -0.545, -0.4, -0.309, -0.109, -0.103, 0.01, 0.119, 0.377,
+        0.79, 0.963, 1.006, 1.115, 1.572, 1.841, 2.047, 2.2
+    ];
+
+    private static readonly IReadOnlyList<double> Y =
+    [
+        80.574, 84.248, 87.264, 87.195, 89.076, 89.608, 89.868, 90.101, 92.405, 95.854, 100.696, 101.06, 401.672,
+        390.724, 567.534, 635.316, 733.054, 759.087, 894.206, 990.785, 1090.109, 1080.914, 1122.643, 1178.351, 1260.531,
+        1273.514, 1288.339, 1327.543, 1353.863, 1414.509, 1425.208, 1421.384, 1442.962, 1464.35, 1468.705, 1447.894,
+        1457.628
+    ];
+
+    private static readonly IReadOnlyList<string> Parameters = ["b1", "b2", "b3", "b4", "b5", "b6", "b7"];
+
+    private static readonly IReadOnlyList<double> InitialValues = [1000, 1000, 400, 40, 0.7, 0.3, 0.03];
+
+    private static readonly IReadOnlyList<double> OptimumValues =
+    [
+        1.2881396800E+03, 1.4910792535E+03, 5.8323836877E+02, 7.5416644291E+01, 9.6629502864E-01, 3.9797285797E-01,
+        4.9727297349E-02
+    ];
+
+    private static readonly Func<double, IReadOnlyList<double>, double> Model = (x, p) =>
+    {
+        var (b1, b2, b3, b4, b5, b6, b7) = (p[0], p[1], p[2], p[3], p[4], p[5], p[6]);
+        return (b1 + b2 * x + b3 * Math.Pow(x, 2) + b4 * Math.Pow(x, 3)) /
+               (b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1);
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelGradient = (x, p) =>
+    {
+        var (b1, b2, b3, b4, b5, b6, b7) = (p[0], p[1], p[2], p[3], p[4], p[5], p[6]);
+        var g0 = 1 / (b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1);
+        var g1 = x / (b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1);
+        var g2 = Math.Pow(x, 2) / (b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1);
+        var g3 = Math.Pow(x, 3) / (b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1);
+        var g4 = -x * (b1 + b2 * x + b3 * Math.Pow(x, 2) + b4 * Math.Pow(x, 3)) /
+                 Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 2);
+        var g5 = Math.Pow(x, 2) * (-b1 - b2 * x - b3 * Math.Pow(x, 2) - b4 * Math.Pow(x, 3)) /
+                 Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 2);
+        var g6 = Math.Pow(x, 3) * (-b1 - b2 * x - b3 * Math.Pow(x, 2) - b4 * Math.Pow(x, 3)) /
+                 Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 2);
+        return [g0, g1, g2, g3, g4, g5, g6];
+    };
+
+    [SuppressMessage("ReSharper", "InlineTemporaryVariable")] 
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessian = (x, p) =>
+    {
+        var (b1, b2, b3, b4, b5, b6, b7) = (p[0], p[1], p[2], p[3], p[4], p[5], p[6]);
+        var h04 = -x / Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 2);
+        var h05 = -Math.Pow(x, 2) / Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 2);
+        var h06 = -Math.Pow(x, 3) / Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 2);
+        var h14 = h05;
+        var h15 = h06;
+        var h16 = -Math.Pow(x, 4) / Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 2);
+        var h24 = h15;
+        var h25 = h16;
+        var h26 = -Math.Pow(x, 5) / Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 2);
+        var h34 = h16;
+        var h35 = h26;
+        var h36 = -Math.Pow(x, 6) / Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 2);
+        var h44 = 2 * Math.Pow(x, 2) * (b1 + b2 * x + b3 * Math.Pow(x, 2) + b4 * Math.Pow(x, 3)) /
+                  Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 3);
+        var h45 = 2 * Math.Pow(x, 3) * (b1 + b2 * x + b3 * Math.Pow(x, 2) + b4 * Math.Pow(x, 3)) /
+                  Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 3);
+        var h46 = 2 * Math.Pow(x, 4) * (b1 + b2 * x + b3 * Math.Pow(x, 2) + b4 * Math.Pow(x, 3)) /
+                  Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 3);
+        var h55 = h46;
+        var h56 = 2 * Math.Pow(x, 5) * (b1 + b2 * x + b3 * Math.Pow(x, 2) + b4 * Math.Pow(x, 3)) /
+                  Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 3);
+        var h66 = 2 * Math.Pow(x, 6) * (b1 + b2 * x + b3 * Math.Pow(x, 2) + b4 * Math.Pow(x, 3)) /
+                  Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 3);
+        return
+        [
+            0.0, 0.0, 0.0, 0.0, h04, h05, h06,
+            0.0, 0.0, 0.0, 0.0, h14, h15, h16,
+            0.0, 0.0, 0.0, 0.0, h24, h25, h26,
+            0.0, 0.0, 0.0, 0.0, h34, h35, h36,
+            h04, h14, h24, h34, h44, h45, h46,
+            h05, h15, h25, h35, h45, h55, h56,
+            h06, h16, h26, h36, h46, h56, h66
+        ];
+    };
+
+    private static readonly Func<double, IReadOnlyList<double>, IReadOnlyList<double>> ModelHessianDiagonal = (x, p) =>
+    {
+        var (b1, b2, b3, b4, b5, b6, b7) = (p[0], p[1], p[2], p[3], p[4], p[5], p[6]);
+        var h44 = 2 * Math.Pow(x, 2) * (b1 + b2 * x + b3 * Math.Pow(x, 2) + b4 * Math.Pow(x, 3)) /
+                  Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 3);
+        var h55 = 2 * Math.Pow(x, 4) * (b1 + b2 * x + b3 * Math.Pow(x, 2) + b4 * Math.Pow(x, 3)) /
+                  Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 3);
+        var h66 = 2 * Math.Pow(x, 6) * (b1 + b2 * x + b3 * Math.Pow(x, 2) + b4 * Math.Pow(x, 3)) /
+                  Math.Pow(b5 * x + b6 * Math.Pow(x, 2) + b7 * Math.Pow(x, 3) + 1, 3);
+        return [0, 0, 0, 0, h44, h55, h66];
+    };
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/RandomExtensions.cs b/test/ExampleProblems/RandomExtensions.cs
new file mode 100644
index 0000000..d821762
--- /dev/null
+++ b/test/ExampleProblems/RandomExtensions.cs
@@ -0,0 +1,14 @@
+namespace ExampleProblems;
+
+internal static class RandomExtensions
+{
+    public static double NextNormal(this Random random, double mean, double standardDeviation)
+    {
+        // Box-Muller transform
+        var u1 = 1.0 - random.NextDouble();
+        var u2 = 1.0 - random.NextDouble();
+        var randomStandardNormal = Math.Sqrt(-2.0 * Math.Log(u1)) * Math.Sin(2.0 * Math.PI * u2);
+
+        return mean + standardDeviation * randomStandardNormal;
+    }
+}
\ No newline at end of file
diff --git a/test/ExampleProblems/Values.cs b/test/ExampleProblems/Values.cs
new file mode 100644
index 0000000..67f4dac
--- /dev/null
+++ b/test/ExampleProblems/Values.cs
@@ -0,0 +1,22 @@
+namespace ExampleProblems;
+
+public static class Values
+{
+    public static double[] LinearlySpacedBetween(double start, double end, double step)
+    {
+        var number = (int)((end - start) / step) + 1;
+        return LinearlySpacedBetween(start, end, number);
+    }
+    
+    private static double[] LinearlySpacedBetween(double start, double end, int number)
+    {
+        var step = (end - start) / (number - 1);
+        return Enumerable.Range(0, number).Select(i => start + i * step).ToArray();
+    }
+
+    public static double[] LogarithmicallySpacedBetween(double start, double end, int number)
+    {
+        var logValues = LinearlySpacedBetween(Math.Log(start), Math.Log(end), number);
+        return logValues.Select(Math.Exp).ToArray();
+    }
+}
\ No newline at end of file
diff --git a/test/minuit2.net.Benchmarks/Benchmarks.cs b/test/minuit2.net.Benchmarks/Benchmarks.cs
deleted file mode 100644
index d1808f9..0000000
--- a/test/minuit2.net.Benchmarks/Benchmarks.cs
+++ /dev/null
@@ -1,86 +0,0 @@
-using BenchmarkDotNet.Attributes;
-using ExampleProblems;
-using minuit2.net.CostFunctions;
-using minuit2.net.Minimizers;
-
-namespace minuit2.net.Benchmarks;
-
-public class Benchmarks
-{
-    private readonly IMinimizer _minimizer = Minimizer.Migrad;
-
-    // By default, problem1 shares all of its parameters with problem2 (since they have the same names);
-    // All other parameters are unique, meaning non-shared, unless renamed
-    private readonly ConfigurableLeastSquaresProblem _problem1 = new QuadraticPolynomialLeastSquaresProblem();
-    private readonly ConfigurableLeastSquaresProblem _problem2 = new CubicPolynomialLeastSquaresProblem();
-    private readonly ConfigurableLeastSquaresProblem _problem3 = new ExponentialDecayLeastSquaresProblem();
-    private readonly ConfigurableLeastSquaresProblem _problem4 = new BellCurveLeastSquaresProblem();
-
-    [Benchmark]
-    public IMinimizationResult BasicMinimizationProblem()
-    {
-        var cost = _problem1.Cost.Build();
-        var parameterConfigurations = _problem1.ParameterConfigurations.Build();
-        return _minimizer.Minimize(cost, parameterConfigurations);
-    }
-
-    [Benchmark]
-    public IMinimizationResult CompositeMinimizationProblemWithoutAnalyticalCostFunctionDerivatives()
-    {
-        var cost = CostFunction.Sum(
-            _problem1.Cost.Build(),
-            _problem2.Cost.Build(),
-            _problem3.Cost.Build(),
-            _problem4.Cost.Build());
-        var parameterConfigurations = _problem2.ParameterConfigurations.Build()
-            .Concat(_problem3.ParameterConfigurations.Build())
-            .Concat(_problem4.ParameterConfigurations.Build())
-            .ToArray();
-        return _minimizer.Minimize(cost, parameterConfigurations);
-    }
-
-    [Benchmark]
-    public IMinimizationResult CompositeMinimizationProblemWithAnalyticalFirstOrderCostFunctionDerivatives()
-    {
-        var cost = CostFunction.Sum(
-            _problem1.Cost.WithGradient().Build(),
-            _problem2.Cost.WithGradient().Build(),
-            _problem3.Cost.WithGradient().Build(),
-            _problem4.Cost.WithGradient().Build());
-        var parameterConfigurations = _problem2.ParameterConfigurations.Build()
-            .Concat(_problem3.ParameterConfigurations.Build())
-            .Concat(_problem4.ParameterConfigurations.Build())
-            .ToArray();
-        return _minimizer.Minimize(cost, parameterConfigurations);
-    }
-
-    [Benchmark]
-    public IMinimizationResult CompositeMinimizationProblemWithAnalyticalSecondOrderCostFunctionDerivatives()
-    {
-        var cost = CostFunction.Sum(
-            _problem1.Cost.WithHessian().Build(),
-            _problem2.Cost.WithHessian().Build(),
-            _problem3.Cost.WithHessian().Build(),
-            _problem4.Cost.WithHessian().Build());
-        var parameterConfigurations = _problem2.ParameterConfigurations.Build()
-            .Concat(_problem3.ParameterConfigurations.Build())
-            .Concat(_problem4.ParameterConfigurations.Build())
-            .ToArray();
-        return _minimizer.Minimize(cost, parameterConfigurations);
-    }
-
-    [Benchmark]
-    public IMinimizationResult CompositeMinimizationProblemWithAnalyticalSecondOrderCostFunctionDerivativesUsingGaussNewtonApproximation()
-    {
-        var cost = CostFunction.Sum(
-            _problem1.Cost.UsingGaussNewtonApproximation().Build(),
-            _problem2.Cost.UsingGaussNewtonApproximation().Build(),
-            _problem3.Cost.UsingGaussNewtonApproximation().Build(),
-            _problem4.Cost.UsingGaussNewtonApproximation().Build());
-        var parameterConfigurations = _problem2.ParameterConfigurations.Build()
-            .Concat(_problem3.ParameterConfigurations.Build())
-            .Concat(_problem4.ParameterConfigurations.Build())
-            .ToArray();
-        return _minimizer.Minimize(cost, parameterConfigurations);
-    }
-}
\ No newline at end of file
diff --git a/test/minuit2.net.Benchmarks/ConfiguredProblemExtensions.cs b/test/minuit2.net.Benchmarks/ConfiguredProblemExtensions.cs
new file mode 100644
index 0000000..88234a6
--- /dev/null
+++ b/test/minuit2.net.Benchmarks/ConfiguredProblemExtensions.cs
@@ -0,0 +1,13 @@
+using ExampleProblems;
+using minuit2.net.Minimizers;
+
+namespace minuit2.net.Benchmarks;
+
+internal static class ConfiguredProblemExtensions
+{
+    public static IMinimizationResult MinimizeWithMigrad(this IConfiguredProblem problem, Strategy strategy)
+    {
+        var minimizerConfiguration = new MinimizerConfiguration(strategy);
+        return Minimizer.Migrad.Minimize(problem.Cost, problem.ParameterConfigurations, minimizerConfiguration);
+    }
+}
\ No newline at end of file
diff --git a/test/minuit2.net.Benchmarks/MinuitTutorialProblemsMigradBenchmarks.cs b/test/minuit2.net.Benchmarks/MinuitTutorialProblemsMigradBenchmarks.cs
new file mode 100644
index 0000000..4cced16
--- /dev/null
+++ b/test/minuit2.net.Benchmarks/MinuitTutorialProblemsMigradBenchmarks.cs
@@ -0,0 +1,39 @@
+using System.Diagnostics.CodeAnalysis;
+using BenchmarkDotNet.Attributes;
+using BenchmarkDotNet.Order;
+using ExampleProblems;
+using ExampleProblems.MinuitTutorialProblems;
+using static ExampleProblems.DerivativeConfiguration;
+
+namespace minuit2.net.Benchmarks;
+
+[SuppressMessage("ReSharper", "UnassignedField.Global", Justification = "[Params] property is set at runtime by Benchmark.NET")]
+[Orderer(SummaryOrderPolicy.Method)]
+public class MinuitTutorialProblemsMigradBenchmarks
+{
+    [Params(WithoutDerivatives, WithGradient, WithGradientAndHessian, WithGradientHessianAndHessianDiagonal)]
+    public DerivativeConfiguration DerivativeConfiguration;
+
+    [Params(Strategy.Fast, Strategy.Balanced, Strategy.Rigorous, Strategy.VeryRigorous)]
+    public Strategy Strategy;
+
+    [Benchmark]
+    public IMinimizationResult RosenbrockProblem() =>
+        new RosenbrockProblem(DerivativeConfiguration).MinimizeWithMigrad(Strategy);
+
+    [Benchmark]
+    public IMinimizationResult WoodProblem() =>
+        new WoodProblem(DerivativeConfiguration).MinimizeWithMigrad(Strategy);
+
+    [Benchmark]
+    public IMinimizationResult PowellProblem() =>
+        new PowellProblem(DerivativeConfiguration).MinimizeWithMigrad(Strategy);
+
+    [Benchmark]
+    public IMinimizationResult FletcherPowellProblem() =>
+        new FletcherPowellProblem(DerivativeConfiguration).MinimizeWithMigrad(Strategy);
+
+    [Benchmark]
+    public IMinimizationResult GoldsteinPriceProblem() =>
+        new GoldsteinPriceProblem(DerivativeConfiguration).MinimizeWithMigrad(Strategy);
+}
\ No newline at end of file
diff --git a/test/minuit2.net.Benchmarks/NistProblemsMigradBenchmarks.cs b/test/minuit2.net.Benchmarks/NistProblemsMigradBenchmarks.cs
new file mode 100644
index 0000000..c52afa9
--- /dev/null
+++ b/test/minuit2.net.Benchmarks/NistProblemsMigradBenchmarks.cs
@@ -0,0 +1,44 @@
+using System.Diagnostics.CodeAnalysis;
+using BenchmarkDotNet.Attributes;
+using BenchmarkDotNet.Order;
+using ExampleProblems;
+using ExampleProblems.NISTProblems;
+using static ExampleProblems.DerivativeConfiguration;
+using static minuit2.net.Strategy;
+
+namespace minuit2.net.Benchmarks;
+
+[SuppressMessage("ReSharper", "UnassignedField.Global", Justification = "[Params] property is set at runtime by Benchmark.NET")]
+[Orderer(SummaryOrderPolicy.Method)]
+public class NistProblemsMigradBenchmarks
+{
+    [Params(WithoutDerivatives, WithGradient, WithGradientAndHessian, WithGradientHessianAndHessianDiagonal)]
+    public DerivativeConfiguration DerivativeConfiguration;
+
+    [Params(Fast, Balanced, Rigorous, VeryRigorous)]
+    public Strategy Strategy;
+
+    [Benchmark]
+    public IMinimizationResult Mgh09Problem() => 
+        new Mgh09Problem(DerivativeConfiguration).MinimizeWithMigrad(Strategy);
+
+    [Benchmark]
+    public IMinimizationResult ThurberProblem() =>
+        new ThurberProblem(DerivativeConfiguration).MinimizeWithMigrad(Strategy);
+    
+    [Benchmark]
+    public IMinimizationResult BoxBodProblem() =>
+        new BoxBodProblem(DerivativeConfiguration).MinimizeWithMigrad(Strategy);
+
+    [Benchmark]
+    public IMinimizationResult Rat42Problem() => 
+        new Rat42Problem(DerivativeConfiguration).MinimizeWithMigrad(Strategy);
+    
+    [Benchmark]
+    public IMinimizationResult Rat43Problem() => 
+        new Rat43Problem(DerivativeConfiguration).MinimizeWithMigrad(Strategy);
+    
+    [Benchmark]
+    public IMinimizationResult Bennett5Problem() => 
+        new Bennett5Problem(DerivativeConfiguration).MinimizeWithMigrad(Strategy);
+}
\ No newline at end of file
diff --git a/test/minuit2.net.Benchmarks/Program.cs b/test/minuit2.net.Benchmarks/Program.cs
index af0ba51..b18f675 100644
--- a/test/minuit2.net.Benchmarks/Program.cs
+++ b/test/minuit2.net.Benchmarks/Program.cs
@@ -1,4 +1,6 @@
 using BenchmarkDotNet.Running;
 using minuit2.net.Benchmarks;
 
-BenchmarkRunner.Run<Benchmarks>();
\ No newline at end of file
+BenchmarkRunner.Run<MinuitTutorialProblemsMigradBenchmarks>();
+BenchmarkRunner.Run<NistProblemsMigradBenchmarks>();
+BenchmarkRunner.Run<SurfaceBiosensorBindingKineticsMigradBenchmarks>();
\ No newline at end of file
diff --git a/test/minuit2.net.Benchmarks/SurfaceBiosensorBindingKineticsMigradBenchmarks.cs b/test/minuit2.net.Benchmarks/SurfaceBiosensorBindingKineticsMigradBenchmarks.cs
new file mode 100644
index 0000000..7952ae9
--- /dev/null
+++ b/test/minuit2.net.Benchmarks/SurfaceBiosensorBindingKineticsMigradBenchmarks.cs
@@ -0,0 +1,34 @@
+using System.Diagnostics.CodeAnalysis;
+using BenchmarkDotNet.Attributes;
+using BenchmarkDotNet.Order;
+using ExampleProblems;
+using ExampleProblems.CustomProblems;
+using static ExampleProblems.DerivativeConfiguration;
+using static minuit2.net.Strategy;
+
+namespace minuit2.net.Benchmarks;
+
+[SuppressMessage("ReSharper", "UnassignedField.Global", Justification = "[Params] property is set at runtime by Benchmark.NET")]
+[Orderer(SummaryOrderPolicy.Method)]
+public class SurfaceBiosensorBindingKineticsMigradBenchmarks
+{
+    [Params(WithoutDerivatives, WithGradient)]
+    // Benchmark all configurations once Hessian-indexing bug is fixed in Minuit2
+    // (see https://github.com/root-project/root/pull/20936)
+    public DerivativeConfiguration DerivativeConfiguration;  
+
+    [Params(Fast, Balanced, Rigorous, VeryRigorous)]
+    public Strategy Strategy;
+
+    [Benchmark]
+    public IMinimizationResult SingleBindingKinetics() =>
+        new SurfaceBiosensorBindingKineticsProblem(DerivativeConfiguration).MinimizeWithMigrad(Strategy);
+
+    [Benchmark]
+    public IMinimizationResult GlobalBindingKinetics()
+    {
+        var analyteConcentrations = Values.LogarithmicallySpacedBetween(1, 100, 50);
+        var problem = SurfaceBiosensorBindingKineticsProblem.Global(analyteConcentrations, DerivativeConfiguration);
+        return problem.MinimizeWithMigrad(Strategy);
+    }
+}
\ No newline at end of file
diff --git a/test/minuit2.net.UnitTests/Any_gradient_based_minimizer.spec.cs b/test/minuit2.net.UnitTests/Any_gradient_based_minimizer.spec.cs
index 9233d48..6c82e75 100644
--- a/test/minuit2.net.UnitTests/Any_gradient_based_minimizer.spec.cs
+++ b/test/minuit2.net.UnitTests/Any_gradient_based_minimizer.spec.cs
@@ -1,6 +1,7 @@
 using AwesomeAssertions;
 using ConstrainedNonDeterminism;
 using ExampleProblems;
+using ExampleProblems.CustomProblems;
 using minuit2.net.CostFunctions;
 using minuit2.net.Minimizers;
 using minuit2.net.UnitTests.TestUtilities;
@@ -11,7 +12,7 @@ namespace minuit2.net.UnitTests;
 public abstract class Any_gradient_based_minimizer(IMinimizer minimizer) : Any_minimizer(minimizer)
 {
     private readonly IMinimizer _minimizer = minimizer;
-    private readonly ConfigurableLeastSquaresProblem _defaultProblem = new CubicPolynomialLeastSquaresProblem();
+    private readonly ConfigurableLeastSquaresProblem _defaultProblem = new CubicPolynomialProblem();
     
     [Test]
     public void when_asked_to_minimize_a_cost_function_with_an_analytical_gradient_that_returns_the_wrong_size_throws_an_exception(
@@ -111,7 +112,7 @@ public void when_minimizing_a_cost_function_with_an_analytical_gradient_yields_a
         var referenceResult = _minimizer.Minimize(referenceCost, parameterConfigurations);
 
         result.Should()
-            .MatchExcludingFunctionCalls(referenceResult, options => options.WithRelativeDoubleTolerance(0.001)).And
+            .MatchExcludingFunctionCalls(referenceResult, options => options.WithSmartDoubleTolerance(0.001)).And
             .HaveFewerFunctionCallsThan(referenceResult);
     }
     
@@ -157,7 +158,7 @@ public void when_minimizing_a_cost_function_with_an_analytical_hessian_yields_a_
         var referenceResult = _minimizer.Minimize(referenceCost, parameterConfigurations);
         
         result.Should()
-            .MatchExcludingFunctionCalls(referenceResult, options => options.WithRelativeDoubleTolerance(0.001)).And
+            .MatchExcludingFunctionCalls(referenceResult, options => options.WithSmartDoubleTolerance(0.001)).And
             .HaveFewerFunctionCallsThan(referenceResult);
     }
     
@@ -201,7 +202,7 @@ public void when_minimizing_a_cost_function_with_an_analytical_hessian_that_is_n
         // Newton step points in the wrong direction — away from the local minimum. To prevent this, the initial
         // Hessian (or its diagonal approximation) must be regularized to ensure positive definiteness during minimizer
         // seeding. Without this safeguard, the minimizer will fail in this case (cf. https://github.com/root-project/root/issues/20665). 
-        var problem = new ExponentialDecayLeastSquaresProblem();
+        var problem = new ExponentialDecayProblem();
         var parameterConfigurations = problem.ParameterConfigurations
             .WithParameter(1).WithLimits(0, null)
             .Build();
@@ -212,7 +213,7 @@ public void when_minimizing_a_cost_function_with_an_analytical_hessian_that_is_n
         var result = _minimizer.Minimize(cost, parameterConfigurations);
         var referenceResult = _minimizer.Minimize(referenceCost, parameterConfigurations);
 
-        result.Should().MatchExcludingFunctionCalls(referenceResult, options => options.WithRelativeDoubleTolerance(0.001));
+        result.Should().MatchExcludingFunctionCalls(referenceResult, options => options.WithSmartDoubleTolerance(0.001));
     }
     
     [Test]
@@ -294,7 +295,7 @@ private IMinimizationResult MinimizeAndRefineErrors(
     public void when_minimizing_a_cost_function_sum_where_only_some_components_have_an_analytical_hessian_yields_the_same_result_as_if_none_had_an_analytical_gradient(
         [Values] Strategy strategy)
     {
-        var problem = new QuadraticPolynomialLeastSquaresProblem();
+        var problem = new QuadraticPolynomialProblem();
         // second component shares offset parameter with the first component
         var component2 = problem.Cost.WithParameterSuffixes("2", [1, 2]).Build();
         var cost = CostFunction.Sum(problem.Cost.WithHessian().Build(), component2);
@@ -315,7 +316,7 @@ public void when_minimizing_a_cost_function_sum_where_all_components_have_analyt
         [Values(1, 2)] double errorDefinitionOfComponent1,
         [Values] Strategy strategy)
     {
-        var problem = new QuadraticPolynomialLeastSquaresProblem();
+        var problem = new QuadraticPolynomialProblem();
         // second component shares offset parameter with the first component
         var cost = CostFunction.Sum(
             problem.Cost.WithErrorDefinition(errorDefinitionOfComponent1).WithHessian().Build(), 
@@ -332,7 +333,7 @@ public void when_minimizing_a_cost_function_sum_where_all_components_have_analyt
         var referenceResult = _minimizer.Minimize(referenceCost, parameterConfigurations, minimizerConfiguration);
         
         result.Should()
-            .MatchExcludingFunctionCalls(referenceResult, options => options.WithRelativeDoubleTolerance(0.001)).And
+            .MatchExcludingFunctionCalls(referenceResult, options => options.WithSmartDoubleTolerance(0.001)).And
             .HaveFewerFunctionCallsThan(referenceResult);
     }
 }
\ No newline at end of file
diff --git a/test/minuit2.net.UnitTests/Any_minimizer.spec.cs b/test/minuit2.net.UnitTests/Any_minimizer.spec.cs
index 6b9a5da..462d1c5 100644
--- a/test/minuit2.net.UnitTests/Any_minimizer.spec.cs
+++ b/test/minuit2.net.UnitTests/Any_minimizer.spec.cs
@@ -1,6 +1,9 @@
 using AwesomeAssertions;
 using ConstrainedNonDeterminism;
 using ExampleProblems;
+using ExampleProblems.CustomProblems;
+using ExampleProblems.MinuitTutorialProblems;
+using ExampleProblems.NISTProblems;
 using minuit2.net.CostFunctions;
 using minuit2.net.Minimizers;
 using minuit2.net.UnitTests.TestUtilities;
@@ -10,21 +13,22 @@ namespace minuit2.net.UnitTests;
 
 public abstract class Any_minimizer(IMinimizer minimizer)
 {
-    private readonly ConfigurableLeastSquaresProblem _defaultProblem = new CubicPolynomialLeastSquaresProblem();
+    private readonly ConfigurableLeastSquaresProblem _defaultProblem = new CubicPolynomialProblem();
     
     protected static IEnumerable<TestCaseData> WellPosedMinimizationProblems()
     {
         foreach (Strategy strategy in Enum.GetValues(typeof(Strategy)))
         {
-            yield return TestCase(new QuadraticPolynomialLeastSquaresProblem().Configured(), nameof(QuadraticPolynomialLeastSquaresProblem));
-            yield return TestCase(new CubicPolynomialLeastSquaresProblem().Configured(), nameof(CubicPolynomialLeastSquaresProblem));
-            yield return TestCase(new ExponentialDecayLeastSquaresProblem().Configured(x => x.WithParameter(1).WithLimits(0, null)), nameof(ExponentialDecayLeastSquaresProblem));
-            yield return TestCase(new BellCurveLeastSquaresProblem().Configured(x => x.WithParameter(1).WithLimits(0, null)), nameof(BellCurveLeastSquaresProblem));
-            yield return TestCase(new NumericalPendulumLeastSquaresProblem(), nameof(NumericalPendulumLeastSquaresProblem));
+            yield return TestCase(new QuadraticPolynomialProblem().Configured(), nameof(QuadraticPolynomialProblem));
+            yield return TestCase(new CubicPolynomialProblem().Configured(), nameof(CubicPolynomialProblem));
+            yield return TestCase(new ExponentialDecayProblem().Configured(x => x.WithParameter(1).WithLimits(0, null)), nameof(ExponentialDecayProblem));
+            yield return TestCase(new BellCurveProblem().Configured(x => x.WithParameter(1).WithLimits(0, null)), nameof(BellCurveProblem));
+            yield return TestCase(new NumericalPendulumProblem(), nameof(NumericalPendulumProblem));
+            yield return TestCase(new SurfaceBiosensorBindingKineticsProblem(), nameof(SurfaceBiosensorBindingKineticsProblem));
             continue;
 
             TestCaseData TestCase(IConfiguredProblem problem, string problemName) =>
-                new TestCaseData(problem, strategy).SetArgDisplayNames(problemName, strategy.ToString());
+                new TestCaseData(problem, strategy).SetName($"{problemName} ({strategy})");
         }
     }
 
@@ -36,8 +40,61 @@ public void when_minimizing_a_well_posed_problem_finds_the_optimum_parameter_val
         var minimizerConfiguration = new MaximumAccuracyMinimizerConfiguration(strategy);
         
         var result = minimizer.Minimize(problem.Cost, problem.ParameterConfigurations, minimizerConfiguration);
+
+        result.ParameterValues.Should().BeApproximately(problem.OptimumParameterValues,
+            relativeToleranceForNonZeros: 0.01, toleranceForZeros: 0.01);
+    }
+    
+    private static IEnumerable<TestCaseData> ChallengingMinimizationProblems()
+    {
+        foreach (Strategy strategy in Enum.GetValues(typeof(Strategy)))
+        foreach (DerivativeConfiguration derivativeConfiguration in Enum.GetValues(typeof(DerivativeConfiguration)))
+        {
+            // Minuit tutorial problems
+            yield return TestCase(new RosenbrockProblem(derivativeConfiguration), nameof(RosenbrockProblem));
+            yield return TestCase(new WoodProblem(derivativeConfiguration), nameof(WoodProblem));
+            yield return TestCase(new PowellProblem(derivativeConfiguration), nameof(PowellProblem));
+            yield return TestCase(new FletcherPowellProblem(derivativeConfiguration), nameof(FletcherPowellProblem));
+            yield return TestCase(new GoldsteinPriceProblem(derivativeConfiguration), nameof(GoldsteinPriceProblem));
+            
+            // NIST problems ranked as difficult 
+            yield return TestCase(new Mgh09Problem(derivativeConfiguration), nameof(Mgh09Problem));
+            yield return TestCase(new ThurberProblem(derivativeConfiguration), nameof(ThurberProblem));
+            yield return TestCase(new BoxBodProblem(derivativeConfiguration), nameof(BoxBodProblem));
+            yield return TestCase(new Rat42Problem(derivativeConfiguration), nameof(Rat42Problem));
+            yield return TestCase(new Rat43Problem(derivativeConfiguration), nameof(Rat43Problem));
+            yield return TestCase(new Bennett5Problem(derivativeConfiguration), nameof(Bennett5Problem));
+            continue;
+
+            TestCaseData TestCase(IConfiguredProblem problem, string problemName) => 
+                new TestCaseData(problem, strategy).SetName($"{problemName}.{derivativeConfiguration} ({strategy})");
+        }
+    }
+
+    [Test]
+    public Task fails_for_the_following_challenging_problems()
+    {
+        var report = new List<string>();
+        foreach (var testCase in ChallengingMinimizationProblems().OrderBy(x => x.TestName))
+        {
+            var problem = (IConfiguredProblem)testCase.Arguments[0]!;
+            var strategy = (Strategy)testCase.Arguments[1]!;
+            
+            var minimizerConfiguration = new MaximumAccuracyMinimizerConfiguration(strategy);
+            var result = minimizer.Minimize(problem.Cost, problem.ParameterConfigurations, minimizerConfiguration);
+
+            try
+            {
+                result.ParameterValues.Should().BeApproximately(problem.OptimumParameterValues, 
+                    relativeToleranceForNonZeros: 0.01, toleranceForZeros: 0.01);
+            }
+            catch (Exception)
+            {
+                report.Add(testCase.TestName!);
+            }
+        }
         
-        result.ParameterValues.Should().BeApproximately(problem.OptimumParameterValues, relativeTolerance: 0.01);
+        return Verify(report);
     }
     
     private static IEnumerable<TestCaseData> InvalidParameterConfigurationTestCases()
@@ -232,7 +289,7 @@ public void when_minimizing_a_cost_function_sum_with_a_single_component_yields_a
         
         sumResult.Should().MatchExcludingFunctionCalls(componentResult, options => options
             .Excluding(x => x.ParameterCovarianceMatrix)  // is null for non-gradient-based minimizers (Simplex)
-            .WithRelativeDoubleTolerance(0.001));
+            .WithSmartDoubleTolerance(0.001));
     }
 
     [Test]
@@ -248,7 +305,7 @@ public void when_minimizing_a_cost_function_sum_of_independent_components_with_d
             [Values] bool hasGradient,
             [Values] Strategy strategy)
     {
-        var problem = new QuadraticPolynomialLeastSquaresProblem();
+        var problem = new QuadraticPolynomialProblem();
         var component1 = problem.Cost.WithParameterSuffixes("1").WithGradient(hasGradient).Build();
         var component2 = problem.Cost.WithParameterSuffixes("2").WithGradient(hasGradient).WithErrorDefinition(2).Build();
         var sum = CostFunction.Sum(component1, component2);
@@ -275,7 +332,7 @@ public void when_minimizing_a_cost_function_sum_of_independent_components_with_d
     public void when_the_cost_function_returns_a_non_finite_value_during_a_minimization_process_yields_an_invalid_result_with_non_finite_value_exit_condition_and_undefined_covariances(
         double nonFiniteValue)
     {
-        var problem = new QuadraticPolynomialLeastSquaresProblem();
+        var problem = new QuadraticPolynomialProblem();
         const int numberOfValidFunctionCalls = 5;
         var cost = problem.Cost.Build().WithValueOverride(_ => nonFiniteValue, numberOfValidFunctionCalls);
         var parameterConfigurations = problem.ParameterConfigurations.Build();
@@ -294,7 +351,7 @@ public void when_the_cost_function_returns_a_non_finite_value_during_a_minimizat
     [Test]
     public void when_the_cost_function_value_calculation_throws_an_exception_during_a_minimization_process_forwards_that_exception()
     {
-        var problem = new QuadraticPolynomialLeastSquaresProblem();
+        var problem = new QuadraticPolynomialProblem();
         var cost = problem.Cost.Build().WithValueOverride(_ => throw new TestException());
         var parameterConfigurations = problem.ParameterConfigurations.Build();
         
diff --git a/test/minuit2.net.UnitTests/Hesse_error_calculator.spec.cs b/test/minuit2.net.UnitTests/Hesse_error_calculator.spec.cs
index 26aa3b7..7040b43 100644
--- a/test/minuit2.net.UnitTests/Hesse_error_calculator.spec.cs
+++ b/test/minuit2.net.UnitTests/Hesse_error_calculator.spec.cs
@@ -1,6 +1,6 @@
 using AwesomeAssertions;
 using ConstrainedNonDeterminism;
-using ExampleProblems;
+using ExampleProblems.CustomProblems;
 using minuit2.net.CostFunctions;
 using minuit2.net.Minimizers;
 using minuit2.net.UnitTests.TestUtilities;
@@ -27,7 +27,7 @@ public class When_refining_a_minimization_result_for_a_related_cost_function
 
         public When_refining_a_minimization_result_for_a_related_cost_function()
         {
-            var problem = new QuadraticPolynomialLeastSquaresProblem();
+            var problem = new QuadraticPolynomialProblem();
             _costFunction = problem.Cost.WithGradient().Build();
             var parameterConfigurations = problem.ParameterConfigurations.Build();
             _minimizationResult = Minimizer.Migrad.Minimize(_costFunction, parameterConfigurations);
@@ -92,7 +92,7 @@ public void When_refining_an_analytical_minimization_result_for_a_related_cost_f
         [Values] bool hasReferenceGradient,
         [Values] Strategy strategy)
     {
-        var problem = new QuadraticPolynomialLeastSquaresProblem();
+        var problem = new QuadraticPolynomialProblem();
         var cost = problem.Cost.WithErrorDefinition(errorDefinition).WithHessian().Build();
         var parameterConfigurations = problem.ParameterConfigurations.Build();
         // minimization result for the (fully) analytical cost function  
@@ -104,7 +104,7 @@ public void When_refining_an_analytical_minimization_result_for_a_related_cost_f
         var referenceResult = HesseErrorCalculator.Refine(analyticalMinimizationResult, referenceCost, strategy);
 
         result.Should()
-            .MatchExcludingFunctionCalls(referenceResult, options => options.WithRelativeDoubleTolerance(0.001)).And
+            .MatchExcludingFunctionCalls(referenceResult, options => options.WithSmartDoubleTolerance(0.001)).And
             .HaveFewerFunctionCallsThan(referenceResult);
     }
 }
\ No newline at end of file
diff --git a/test/minuit2.net.UnitTests/Least_squares_cost_function.spec.cs b/test/minuit2.net.UnitTests/Least_squares_cost_function.spec.cs
index 0baf972..7ba6b0c 100644
--- a/test/minuit2.net.UnitTests/Least_squares_cost_function.spec.cs
+++ b/test/minuit2.net.UnitTests/Least_squares_cost_function.spec.cs
@@ -106,7 +106,7 @@ public void and_an_analytical_model_gradient_when_asked_for_its_gradient_returns
             }
 
             cost.GradientFor(_parameterValues).Should().BeEquivalentTo(expectedGradient, 
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
         
         [Test]
@@ -126,7 +126,7 @@ public void and_analytical_model_gradient_and_hessian_when_asked_for_its_hessian
             }
 
             cost.HessianFor(_parameterValues).Should().BeEquivalentTo(expectedHessian, 
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
         
         [Test]
@@ -144,7 +144,7 @@ public void and_analytical_model_gradient_and_hessian_when_asked_for_its_hessian
             }
 
             cost.HessianDiagonalFor(_parameterValues).Should().BeEquivalentTo(expectedHessianDiagonal, 
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
         
         [Test]
@@ -162,7 +162,7 @@ public void and_analytical_model_gradient_and_hessian_and_hessian_diagonal_when_
             }
 
             cost.HessianDiagonalFor(_parameterValues).Should().BeEquivalentTo(expectedHessianDiagonal, 
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
         
         [Test, 
@@ -286,7 +286,7 @@ public void and_an_analytical_model_gradient_when_asked_for_its_gradient_returns
             }
             
             cost.GradientFor(_parameterValues).Should().BeEquivalentTo(expectedGradient, 
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
         
         [Test]
@@ -306,7 +306,7 @@ public void and_analytical_model_gradient_and_hessian_when_asked_for_its_hessian
             }
 
             cost.HessianFor(_parameterValues).Should().BeEquivalentTo(expectedHessian, 
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
         
         [Test]
@@ -324,7 +324,7 @@ public void and_analytical_model_gradient_and_hessian_when_asked_for_its_hessian
             }
 
             cost.HessianDiagonalFor(_parameterValues).Should().BeEquivalentTo(expectedHessianDiagonal, 
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
         
         [Test]
@@ -342,7 +342,7 @@ public void and_analytical_model_gradient_and_hessian_and_hessian_diagonal_when_
             }
 
             cost.HessianDiagonalFor(_parameterValues).Should().BeEquivalentTo(expectedHessianDiagonal, 
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
         
         [Test, 
@@ -454,7 +454,7 @@ public void and_an_analytical_model_gradient_when_asked_for_its_gradient_returns
             }
             
             cost.GradientFor(_parameterValues).Should().BeEquivalentTo(expectedGradient, 
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
         
         [Test]
@@ -474,7 +474,7 @@ public void and_analytical_model_gradient_and_hessian_when_asked_for_its_hessian
             }
 
             cost.HessianFor(_parameterValues).Should().BeEquivalentTo(expectedHessian, 
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
         
         [Test]
@@ -492,7 +492,7 @@ public void and_analytical_model_gradient_and_hessian_when_asked_for_its_hessian
             }
 
             cost.HessianDiagonalFor(_parameterValues).Should().BeEquivalentTo(expectedHessianDiagonal, 
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
         
         [Test]
@@ -510,7 +510,7 @@ public void and_analytical_model_gradient_and_hessian_and_hessian_diagonal_when_
             }
 
             cost.HessianDiagonalFor(_parameterValues).Should().BeEquivalentTo(expectedHessianDiagonal, 
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
         
         [Test, 
diff --git a/test/minuit2.net.UnitTests/Least_squares_cost_function_with_gauss_newton_approximation.spec.cs b/test/minuit2.net.UnitTests/Least_squares_cost_function_with_gauss_newton_approximation.spec.cs
index b026062..cde7a94 100644
--- a/test/minuit2.net.UnitTests/Least_squares_cost_function_with_gauss_newton_approximation.spec.cs
+++ b/test/minuit2.net.UnitTests/Least_squares_cost_function_with_gauss_newton_approximation.spec.cs
@@ -110,7 +110,7 @@ public void when_asked_for_its_gradient_returns_the_expected_vector()
             }
 
             _validCost.GradientFor(_parameterValues).Should().BeEquivalentTo(expectedGradient,
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
 
         [Test]
@@ -127,7 +127,7 @@ public void when_asked_for_its_hessian_returns_the_expected_flat_matrix()
             }
 
             _validCost.HessianFor(_parameterValues).Should().BeEquivalentTo(expectedHessian,
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
 
         [Test]
@@ -142,7 +142,7 @@ public void when_asked_for_its_hessian_diagonal_returns_the_expected_vector()
             }
 
             _validCost.HessianDiagonalFor(_parameterValues).Should().BeEquivalentTo(expectedHessianDiagonal,
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
     }
 
@@ -254,7 +254,7 @@ public void when_asked_for_its_gradient_returns_the_expected_vector()
             }
 
             _validCost.GradientFor(_parameterValues).Should().BeEquivalentTo(expectedGradient,
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
 
         [Test]
@@ -271,7 +271,7 @@ public void when_asked_for_its_hessian_returns_the_expected_flat_matrix()
             }
 
             _validCost.HessianFor(_parameterValues).Should().BeEquivalentTo(expectedHessian,
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
 
         [Test]
@@ -286,7 +286,7 @@ public void when_asked_for_its_hessian_diagonal_returns_the_expected_vector()
             }
 
             _validCost.HessianDiagonalFor(_parameterValues).Should().BeEquivalentTo(expectedHessianDiagonal,
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
     }
 
@@ -385,7 +385,7 @@ public void when_asked_for_its_gradient_returns_the_expected_vector()
             }
 
             _validCost.GradientFor(_parameterValues).Should().BeEquivalentTo(expectedGradient,
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
 
         [Test]
@@ -402,7 +402,7 @@ public void when_asked_for_its_hessian_returns_the_expected_flat_matrix()
             }
 
             _validCost.HessianFor(_parameterValues).Should().BeEquivalentTo(expectedHessian,
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
 
         [Test]
@@ -417,7 +417,7 @@ public void when_asked_for_its_hessian_diagonal_returns_the_expected_vector()
             }
 
             _validCost.HessianDiagonalFor(_parameterValues).Should().BeEquivalentTo(expectedHessianDiagonal,
-                options => options.WithRelativeDoubleTolerance(0.001));
+                options => options.WithSmartDoubleTolerance(0.001));
         }
     }
 }
\ No newline at end of file
diff --git a/test/minuit2.net.UnitTests/Migrad_minimizer.spec.cs b/test/minuit2.net.UnitTests/Migrad_minimizer.spec.cs
index 703e2f2..de530f1 100644
--- a/test/minuit2.net.UnitTests/Migrad_minimizer.spec.cs
+++ b/test/minuit2.net.UnitTests/Migrad_minimizer.spec.cs
@@ -1,5 +1,6 @@
 using AwesomeAssertions;
 using ExampleProblems;
+using ExampleProblems.CustomProblems;
 using minuit2.net.CostFunctions;
 using minuit2.net.Minimizers;
 using minuit2.net.UnitTests.TestUtilities;
@@ -25,7 +26,7 @@ private static IMinimizationResult MinimizeAndRefineErrors(
     
     public class When_minimizing_a_least_squares_cost_function
     {
-        private readonly ConfigurableLeastSquaresProblem _problem = new CubicPolynomialLeastSquaresProblem();
+        private readonly ConfigurableLeastSquaresProblem _problem = new CubicPolynomialProblem();
 
         [TestCase(false, 100)]
         [TestCase(true, 78)]
@@ -213,7 +214,7 @@ public void with_unknown_data_errors_and_subsequently_applying_error_definition_
                      "state with clear error details.")]
         public void leaving_the_valid_parameter_space_during_the_minimization_process_yields_an_invalid_result_with_non_finite_value_exit_condition_and_undefined_covariances_representing_the_last_state_of_the_process()
         {
-            var problem = new BellCurveLeastSquaresProblem();
+            var problem = new BellCurveProblem();
             var cost = problem.Cost.Build();
             var parameterConfigurations = problem.ParameterConfigurations
                 .WithParameter(0).WithValue(6.41).And
@@ -235,7 +236,7 @@ public void leaving_the_valid_parameter_space_during_the_minimization_process_yi
 
     public class When_minimizing_a_sum_of_least_squares_cost_functions
     {
-        private readonly ConfigurableLeastSquaresProblem _problem = new CubicPolynomialLeastSquaresProblem();
+        private readonly ConfigurableLeastSquaresProblem _problem = new CubicPolynomialProblem();
 
         [Test]
         [Description("Ensures that recalculation of the error definition for least squares cost functions with " +
diff --git a/test/minuit2.net.UnitTests/TestUtilities/EquivalencyOptionsExtensions.cs b/test/minuit2.net.UnitTests/TestUtilities/EquivalencyOptionsExtensions.cs
index 7bb8b7c..280d457 100644
--- a/test/minuit2.net.UnitTests/TestUtilities/EquivalencyOptionsExtensions.cs
+++ b/test/minuit2.net.UnitTests/TestUtilities/EquivalencyOptionsExtensions.cs
@@ -6,29 +6,64 @@ namespace minuit2.net.UnitTests.TestUtilities;
 
 internal static class EquivalencyOptionsExtensions
 {
-    public static EquivalencyOptions<double> WithRelativeDoubleTolerance(
+    public static EquivalencyOptions<double> WithSmartDoubleTolerance(
         this EquivalencyOptions<double> options,
-        double relativeTolerance)
+        double relativeToleranceForNonZeros, 
+        double? minimumToleranceForNonZeros = null,
+        double? toleranceForZeros = null)
     {
-        return options.WithRelativeDoubleTolerance<double>(relativeTolerance);
+        return options.WithSmartDoubleTolerance<double>(
+            relativeToleranceForNonZeros,
+            minimumToleranceForNonZeros,
+            toleranceForZeros);
     }
     
-    public static EquivalencyOptions<double[,]> WithRelativeDoubleTolerance(
-        this EquivalencyOptions<double[,]> options, double relativeTolerance)
+    public static EquivalencyOptions<double[,]> WithSmartDoubleTolerance(
+        this EquivalencyOptions<double[,]> options, 
+        double relativeToleranceForNonZeros,
+        double? minimumToleranceForNonZeros = null,
+        double? toleranceForZeros = null)
     {
-        return options.WithRelativeDoubleTolerance<double[,]>(relativeTolerance);
+        return options.WithSmartDoubleTolerance<double[,]>(
+            relativeToleranceForNonZeros,
+            minimumToleranceForNonZeros,
+            toleranceForZeros);
     }
     
-    public static EquivalencyOptions<IMinimizationResult> WithRelativeDoubleTolerance(
-        this EquivalencyOptions<IMinimizationResult> options, double relativeTolerance)
+    public static EquivalencyOptions<IMinimizationResult> WithSmartDoubleTolerance(
+        this EquivalencyOptions<IMinimizationResult> options, 
+        double relativeToleranceForNonZeros, 
+        double? minimumToleranceForNonZeros = null,
+        double? toleranceForZeros = null)
     {
-        return options.WithRelativeDoubleTolerance<IMinimizationResult>(relativeTolerance);
+        return options.WithSmartDoubleTolerance<IMinimizationResult>(
+            relativeToleranceForNonZeros,
+            minimumToleranceForNonZeros,
+            toleranceForZeros);
     }
 
-    private static EquivalencyOptions<T> WithRelativeDoubleTolerance<T>(this EquivalencyOptions<T> options, double relativeTolerance)
+    private static EquivalencyOptions<T> WithSmartDoubleTolerance<T>(
+        this EquivalencyOptions<T> options,
+        double relativeToleranceForNonZeros,
+        double? minimumToleranceForNonZeros,
+        double? toleranceForZeros)
     {
         return options
-            .Using<double>(ctx => ctx.Subject.Should().BeApproximately(ctx.Expectation, relativeTolerance, DefaultMinimumDoubleTolerance))
+            .Using<double>(ctx =>
+            {
+                var actual = ctx.Subject;
+                var expected = ctx.Expectation;
+                if (Math.Abs(expected) > double.Epsilon)
+                {
+                    var minimumTolerance = minimumToleranceForNonZeros ?? DefaultDoubleTolerance;
+                    actual.Should().BeApproximately(expected, relativeToleranceForNonZeros, minimumTolerance);
+                }
+                else
+                {
+                    var tolerance = toleranceForZeros ?? DefaultDoubleTolerance;
+                    actual.Should().BeApproximately(expected, tolerance);
+                }
+            })
             .WhenTypeIs<double>();
     }
 }
\ No newline at end of file
diff --git a/test/minuit2.net.UnitTests/TestUtilities/GenericCollectionAssertionExtensions.cs b/test/minuit2.net.UnitTests/TestUtilities/GenericCollectionAssertionExtensions.cs
index 2a71066..65d89f3 100644
--- a/test/minuit2.net.UnitTests/TestUtilities/GenericCollectionAssertionExtensions.cs
+++ b/test/minuit2.net.UnitTests/TestUtilities/GenericCollectionAssertionExtensions.cs
@@ -14,19 +14,27 @@ public static AndConstraint<GenericCollectionAssertions<double>> BeApproximately
         [StringSyntax("CompositeFormat")] string because = "",
         params object[] becauseArgs)
     {
-        return parent.BeApproximately(expectation, DefaultRelativeDoubleTolerance, because, becauseArgs);
+        return parent.BeApproximately(
+            expectation,
+            DefaultRelativeDoubleTolerance,
+            DefaultDoubleTolerance,
+            DefaultDoubleTolerance,
+            because,
+            becauseArgs);
     }
     
     public static AndConstraint<GenericCollectionAssertions<double>> BeApproximately(
         this GenericCollectionAssertions<double> parent,
         IEnumerable<double> expectation,
-        double relativeTolerance,
+        double relativeToleranceForNonZeros,
+        double? minimumToleranceForNonZeros = null,
+        double? toleranceForZeros = null,
         [StringSyntax("CompositeFormat")] string because = "",
         params object[] becauseArgs)
     {
         return parent.BeEquivalentTo(
             expectation, 
-            options => options.WithRelativeDoubleTolerance(relativeTolerance), 
+            options => options.WithSmartDoubleTolerance(relativeToleranceForNonZeros, minimumToleranceForNonZeros, toleranceForZeros), 
             because, 
             becauseArgs);
     }
diff --git a/test/minuit2.net.UnitTests/TestUtilities/NumericAssertionExtensions.cs b/test/minuit2.net.UnitTests/TestUtilities/NumericAssertionExtensions.cs
index 833059f..5553858 100644
--- a/test/minuit2.net.UnitTests/TestUtilities/NumericAssertionExtensions.cs
+++ b/test/minuit2.net.UnitTests/TestUtilities/NumericAssertionExtensions.cs
@@ -8,7 +8,7 @@ namespace minuit2.net.UnitTests.TestUtilities;
 internal static class NumericAssertionExtensions
 {
     internal const double DefaultRelativeDoubleTolerance = 0.001;
-    internal const double DefaultMinimumDoubleTolerance = 1E-8;
+    internal const double DefaultDoubleTolerance = 1E-8;
     
     public static AndConstraint<NumericAssertions<double>> BeApproximately(
         this NumericAssertions<double> parent, 
@@ -19,7 +19,7 @@ public static AndConstraint<NumericAssertions<double>> BeApproximately(
         return parent.BeApproximately(
             expectedValue, 
             DefaultRelativeDoubleTolerance, 
-            DefaultMinimumDoubleTolerance, 
+            DefaultDoubleTolerance, 
             because, 
             becauseArgs);
     }
diff --git a/test/minuit2.net.UnitTests/TestUtilities/ObjectAssertionExtensions.cs b/test/minuit2.net.UnitTests/TestUtilities/ObjectAssertionExtensions.cs
index 6998cfa..6ae33d9 100644
--- a/test/minuit2.net.UnitTests/TestUtilities/ObjectAssertionExtensions.cs
+++ b/test/minuit2.net.UnitTests/TestUtilities/ObjectAssertionExtensions.cs
@@ -31,7 +31,7 @@ public static AndConstraint<ObjectAssertions> BeApproximately(
             ? parent.BeNull()
             : parent.BeEquivalentTo(
                 expectation,
-                options => options.WithRelativeDoubleTolerance(relativeTolerance),
+                options => options.WithSmartDoubleTolerance(relativeTolerance),
                 because,
                 becauseArgs);
     }
diff --git a/test/minuit2.net.UnitTests/The_combined_minimizer.fails_for_the_following_challenging_problems.verified.txt b/test/minuit2.net.UnitTests/The_combined_minimizer.fails_for_the_following_challenging_problems.verified.txt
new file mode 100644
index 0000000..f3b0453
--- /dev/null
+++ b/test/minuit2.net.UnitTests/The_combined_minimizer.fails_for_the_following_challenging_problems.verified.txt
@@ -0,0 +1,35 @@
+[
+  Bennett5Problem.WithGradient (Balanced),
+  Bennett5Problem.WithGradient (Fast),
+  Bennett5Problem.WithGradient (Rigorous),
+  Bennett5Problem.WithGradient (VeryRigorous),
+  Bennett5Problem.WithoutDerivatives (Balanced),
+  Bennett5Problem.WithoutDerivatives (Fast),
+  Bennett5Problem.WithoutDerivatives (Rigorous),
+  Bennett5Problem.WithoutDerivatives (VeryRigorous),
+  BoxBodProblem.WithGradient (Rigorous),
+  BoxBodProblem.WithGradientAndHessian (Rigorous),
+  BoxBodProblem.WithGradientHessianAndHessianDiagonal (Rigorous),
+  BoxBodProblem.WithoutDerivatives (Rigorous),
+  GoldsteinPriceProblem.WithGradient (Balanced),
+  GoldsteinPriceProblem.WithGradient (Fast),
+  GoldsteinPriceProblem.WithGradient (Rigorous),
+  GoldsteinPriceProblem.WithGradient (VeryRigorous),
+  GoldsteinPriceProblem.WithGradientAndHessian (Balanced),
+  GoldsteinPriceProblem.WithGradientAndHessian (Fast),
+  GoldsteinPriceProblem.WithGradientAndHessian (Rigorous),
+  GoldsteinPriceProblem.WithGradientAndHessian (VeryRigorous),
+  GoldsteinPriceProblem.WithGradientHessianAndHessianDiagonal (Balanced),
+  GoldsteinPriceProblem.WithGradientHessianAndHessianDiagonal (Fast),
+  GoldsteinPriceProblem.WithGradientHessianAndHessianDiagonal (Rigorous),
+  GoldsteinPriceProblem.WithGradientHessianAndHessianDiagonal (VeryRigorous),
+  GoldsteinPriceProblem.WithoutDerivatives (Balanced),
+  GoldsteinPriceProblem.WithoutDerivatives (Fast),
+  GoldsteinPriceProblem.WithoutDerivatives (Rigorous),
+  GoldsteinPriceProblem.WithoutDerivatives (VeryRigorous),
+  Mgh09Problem.WithGradient (Balanced),
+  Mgh09Problem.WithGradient (Fast),
+  Mgh09Problem.WithGradientAndHessian (Fast),
+  Mgh09Problem.WithGradientHessianAndHessianDiagonal (Fast),
+  Mgh09Problem.WithoutDerivatives (Fast)
+]
\ No newline at end of file
diff --git a/test/minuit2.net.UnitTests/The_migrad_minimizer.fails_for_the_following_challenging_problems.verified.txt b/test/minuit2.net.UnitTests/The_migrad_minimizer.fails_for_the_following_challenging_problems.verified.txt
new file mode 100644
index 0000000..18d06b9
--- /dev/null
+++ b/test/minuit2.net.UnitTests/The_migrad_minimizer.fails_for_the_following_challenging_problems.verified.txt
@@ -0,0 +1,37 @@
+[
+  Bennett5Problem.WithGradient (Balanced),
+  Bennett5Problem.WithGradient (Fast),
+  Bennett5Problem.WithGradient (Rigorous),
+  Bennett5Problem.WithGradient (VeryRigorous),
+  Bennett5Problem.WithoutDerivatives (Balanced),
+  Bennett5Problem.WithoutDerivatives (Fast),
+  Bennett5Problem.WithoutDerivatives (Rigorous),
+  Bennett5Problem.WithoutDerivatives (VeryRigorous),
+  BoxBodProblem.WithGradient (Rigorous),
+  BoxBodProblem.WithGradientAndHessian (Rigorous),
+  BoxBodProblem.WithGradientHessianAndHessianDiagonal (Rigorous),
+  BoxBodProblem.WithoutDerivatives (Rigorous),
+  GoldsteinPriceProblem.WithGradient (Balanced),
+  GoldsteinPriceProblem.WithGradient (Fast),
+  GoldsteinPriceProblem.WithGradient (Rigorous),
+  GoldsteinPriceProblem.WithGradient (VeryRigorous),
+  GoldsteinPriceProblem.WithGradientAndHessian (Balanced),
+  GoldsteinPriceProblem.WithGradientAndHessian (Fast),
+  GoldsteinPriceProblem.WithGradientAndHessian (Rigorous),
+  GoldsteinPriceProblem.WithGradientAndHessian (VeryRigorous),
+  GoldsteinPriceProblem.WithGradientHessianAndHessianDiagonal (Balanced),
+  GoldsteinPriceProblem.WithGradientHessianAndHessianDiagonal (Fast),
+  GoldsteinPriceProblem.WithGradientHessianAndHessianDiagonal (Rigorous),
+  GoldsteinPriceProblem.WithGradientHessianAndHessianDiagonal (VeryRigorous),
+  GoldsteinPriceProblem.WithoutDerivatives (Balanced),
+  GoldsteinPriceProblem.WithoutDerivatives (Fast),
+  GoldsteinPriceProblem.WithoutDerivatives (Rigorous),
+  GoldsteinPriceProblem.WithoutDerivatives (VeryRigorous),
+  Mgh09Problem.WithGradient (Balanced),
+  Mgh09Problem.WithGradient (Fast),
+  Mgh09Problem.WithGradientAndHessian (Fast),
+  Mgh09Problem.WithGradientHessianAndHessianDiagonal (Fast),
+  Mgh09Problem.WithoutDerivatives (Balanced),
+  Mgh09Problem.WithoutDerivatives (Fast),
+  Mgh09Problem.WithoutDerivatives (Rigorous)
+]
\ No newline at end of file
diff --git a/test/minuit2.net.UnitTests/The_simplex_minimizer.fails_for_the_following_challenging_problems.verified.txt b/test/minuit2.net.UnitTests/The_simplex_minimizer.fails_for_the_following_challenging_problems.verified.txt
new file mode 100644
index 0000000..66855cf
--- /dev/null
+++ b/test/minuit2.net.UnitTests/The_simplex_minimizer.fails_for_the_following_challenging_problems.verified.txt
@@ -0,0 +1,98 @@
+[
+  Bennett5Problem.WithGradient (Balanced),
+  Bennett5Problem.WithGradient (Fast),
+  Bennett5Problem.WithGradient (Rigorous),
+  Bennett5Problem.WithGradient (VeryRigorous),
+  Bennett5Problem.WithGradientAndHessian (Balanced),
+  Bennett5Problem.WithGradientAndHessian (Fast),
+  Bennett5Problem.WithGradientAndHessian (Rigorous),
+  Bennett5Problem.WithGradientAndHessian (VeryRigorous),
+  Bennett5Problem.WithGradientHessianAndHessianDiagonal (Balanced),
+  Bennett5Problem.WithGradientHessianAndHessianDiagonal (Fast),
+  Bennett5Problem.WithGradientHessianAndHessianDiagonal (Rigorous),
+  Bennett5Problem.WithGradientHessianAndHessianDiagonal (VeryRigorous),
+  Bennett5Problem.WithoutDerivatives (Balanced),
+  Bennett5Problem.WithoutDerivatives (Fast),
+  Bennett5Problem.WithoutDerivatives (Rigorous),
+  Bennett5Problem.WithoutDerivatives (VeryRigorous),
+  BoxBodProblem.WithGradient (Balanced),
+  BoxBodProblem.WithGradient (Fast),
+  BoxBodProblem.WithGradient (Rigorous),
+  BoxBodProblem.WithGradient (VeryRigorous),
+  BoxBodProblem.WithGradientAndHessian (Balanced),
+  BoxBodProblem.WithGradientAndHessian (Fast),
+  BoxBodProblem.WithGradientAndHessian (Rigorous),
+  BoxBodProblem.WithGradientAndHessian (VeryRigorous),
+  BoxBodProblem.WithGradientHessianAndHessianDiagonal (Balanced),
+  BoxBodProblem.WithGradientHessianAndHessianDiagonal (Fast),
+  BoxBodProblem.WithGradientHessianAndHessianDiagonal (Rigorous),
+  BoxBodProblem.WithGradientHessianAndHessianDiagonal (VeryRigorous),
+  BoxBodProblem.WithoutDerivatives (Balanced),
+  BoxBodProblem.WithoutDerivatives (Fast),
+  BoxBodProblem.WithoutDerivatives (Rigorous),
+  BoxBodProblem.WithoutDerivatives (VeryRigorous),
+  FletcherPowellProblem.WithGradient (Balanced),
+  FletcherPowellProblem.WithGradient (Fast),
+  FletcherPowellProblem.WithGradient (Rigorous),
+  FletcherPowellProblem.WithGradient (VeryRigorous),
+  FletcherPowellProblem.WithGradientAndHessian (Balanced),
+  FletcherPowellProblem.WithGradientAndHessian (Fast),
+  FletcherPowellProblem.WithGradientAndHessian (Rigorous),
+  FletcherPowellProblem.WithGradientAndHessian (VeryRigorous),
+  FletcherPowellProblem.WithGradientHessianAndHessianDiagonal (Balanced),
+  FletcherPowellProblem.WithGradientHessianAndHessianDiagonal (Fast),
+  FletcherPowellProblem.WithGradientHessianAndHessianDiagonal (Rigorous),
+  FletcherPowellProblem.WithGradientHessianAndHessianDiagonal (VeryRigorous),
+  FletcherPowellProblem.WithoutDerivatives (Balanced),
+  FletcherPowellProblem.WithoutDerivatives (Fast),
+  FletcherPowellProblem.WithoutDerivatives (Rigorous),
+  FletcherPowellProblem.WithoutDerivatives (VeryRigorous),
+  Mgh09Problem.WithGradient (Balanced),
+  Mgh09Problem.WithGradient (Fast),
+  Mgh09Problem.WithGradient (Rigorous),
+  Mgh09Problem.WithGradient (VeryRigorous),
+  Mgh09Problem.WithGradientAndHessian (Balanced),
+  Mgh09Problem.WithGradientAndHessian (Fast),
+  Mgh09Problem.WithGradientAndHessian (Rigorous),
+  Mgh09Problem.WithGradientAndHessian (VeryRigorous),
+  Mgh09Problem.WithGradientHessianAndHessianDiagonal (Balanced),
+  Mgh09Problem.WithGradientHessianAndHessianDiagonal (Fast),
+  Mgh09Problem.WithGradientHessianAndHessianDiagonal (Rigorous),
+  Mgh09Problem.WithGradientHessianAndHessianDiagonal (VeryRigorous),
+  Mgh09Problem.WithoutDerivatives (Balanced),
+  Mgh09Problem.WithoutDerivatives (Fast),
+  Mgh09Problem.WithoutDerivatives (Rigorous),
+  Mgh09Problem.WithoutDerivatives (VeryRigorous),
+  Rat43Problem.WithGradient (Balanced),
+  Rat43Problem.WithGradient (Fast),
+  Rat43Problem.WithGradient (Rigorous),
+  Rat43Problem.WithGradient (VeryRigorous),
+  Rat43Problem.WithGradientAndHessian (Balanced),
+  Rat43Problem.WithGradientAndHessian (Fast),
+  Rat43Problem.WithGradientAndHessian (Rigorous),
+  Rat43Problem.WithGradientAndHessian (VeryRigorous),
+  Rat43Problem.WithGradientHessianAndHessianDiagonal (Balanced),
+  Rat43Problem.WithGradientHessianAndHessianDiagonal (Fast),
+  Rat43Problem.WithGradientHessianAndHessianDiagonal (Rigorous),
+  Rat43Problem.WithGradientHessianAndHessianDiagonal (VeryRigorous),
+  Rat43Problem.WithoutDerivatives (Balanced),
+  Rat43Problem.WithoutDerivatives (Fast),
+  Rat43Problem.WithoutDerivatives (Rigorous),
+  Rat43Problem.WithoutDerivatives (VeryRigorous),
+  ThurberProblem.WithGradient (Balanced),
+  ThurberProblem.WithGradient (Fast),
+  ThurberProblem.WithGradient (Rigorous),
+  ThurberProblem.WithGradient (VeryRigorous),
+  ThurberProblem.WithGradientAndHessian (Balanced),
+  ThurberProblem.WithGradientAndHessian (Fast),
+  ThurberProblem.WithGradientAndHessian (Rigorous),
+  ThurberProblem.WithGradientAndHessian (VeryRigorous),
+  ThurberProblem.WithGradientHessianAndHessianDiagonal (Balanced),
+  ThurberProblem.WithGradientHessianAndHessianDiagonal (Fast),
+  ThurberProblem.WithGradientHessianAndHessianDiagonal (Rigorous),
+  ThurberProblem.WithGradientHessianAndHessianDiagonal (VeryRigorous),
+  ThurberProblem.WithoutDerivatives (Balanced),
+  ThurberProblem.WithoutDerivatives (Fast),
+  ThurberProblem.WithoutDerivatives (Rigorous),
+  ThurberProblem.WithoutDerivatives (VeryRigorous)
+]
\ No newline at end of file
diff --git a/test/minuit2.net.UnitTests/minuit2.net.UnitTests.csproj b/test/minuit2.net.UnitTests/minuit2.net.UnitTests.csproj
index eb41f1c..3dd2a09 100644
--- a/test/minuit2.net.UnitTests/minuit2.net.UnitTests.csproj
+++ b/test/minuit2.net.UnitTests/minuit2.net.UnitTests.csproj
@@ -19,9 +19,10 @@
         <PackageReference Include="AwesomeAssertions" Version="9.1.0" />
         <PackageReference Include="coverlet.collector" Version="6.0.0"/>
         <PackageReference Include="Microsoft.NET.Test.Sdk" Version="17.8.0"/>
-        <PackageReference Include="NUnit" Version="3.14.0"/>
-        <PackageReference Include="NUnit.Analyzers" Version="3.9.0"/>
-        <PackageReference Include="NUnit3TestAdapter" Version="4.5.0"/>
+        <PackageReference Include="NUnit" Version="4.4.0" />
+        <PackageReference Include="NUnit.Analyzers" Version="4.11.2"/>
+        <PackageReference Include="NUnit3TestAdapter" Version="6.1.0"/>
+        <PackageReference Include="Verify.NUnit" Version="31.9.4" />
     </ItemGroup>
 
     <ItemGroup>