Почему не работает метод кубической параболы для поиска минимума функции на интервале? C++

Данный код на отрезке [4,8] должен вывести минимум 28.4565, а выводится -4.00505:

#include <iostream>
#include <cmath>

using namespace std;

double f1(double);
double d1(double);

double cubicParabola(double, double, double, double (*)(double), double (*)(double), int &);

int main(){

double a, b, h, x0;
 double epsilon = pow(10, -6);
 cout << " Функция: x^2 - 4x*sin(x) + cos(x)" << endl;
  cout << "Введите левую границу интервала: " << endl;
  cin >> a;
  cout << "Введите правую границу интервала: " << endl;
  cin >> b;
 double minXCubic = cubicParabola(a, b, epsilon, *f, *d);
 double minCubic = f(minXCubic);
 cout << "Точка минимума: " << minXCubic << endl;
 cout << "Минимум: " << minCubic << endl;
 cout << "Количество итераций: " << iterCount2 << endl;

return 0;
}

double f1(double x) {
    return x * x - 4 * x * sin(x) + cos(x); 
}

double d1(double x) {
    return 2 * x - 5 * sin(x) - 4 * x * cos(x);
}

double cubicParabola(double a, double b, double epsilon, double (*f)(double), double (*d)(double))
{
    int m1 = 100;
    double x1, x2, x3, h1, z1;
    double zm = 0;
    double D1, D2, Y1, Y2, Y3, r, p, q, podkor;
    h1 = (b - a) / m1;
    x1 = a;
    D1 = d(x1);
    if (D1 < 0)
        h1 = -h1;
    x2 = x1 + h1;
    D2 = d(x2);
    Y1 =  f(x1);
    Y2 = f(x2);
   if ((D2 - D1) / h1 < 0) {
        do
        {
            z1 = x1 - x2;
            r = D2;
            p = (D1 - D2 - 2 * (Y1 - Y2 - D2 * z1) / z1) / (z1 * z1);
            q = (D2 - D1 + 3 * (Y1 - Y2 - D2 * z1) / z1) / z1;
            podkor = q * q - 3 * p * r;
            if (podkor < 0)
                zm = -q / (3 * p);
            else
                zm = (-q + sqrt(podkor)) / (3 * p);
            x1 = x2;
            Y1 = Y2;
            D1 = D2;
            x2 = x2 + zm;
            Y2 = f(x2);
            D2 = d(x2);
        } while (abs(zm) > epsilon);
    }
   else {
       cout << "Начальное приближение выбрано неправильно" << endl;
       return 0;
   }
    x3 = x2 + zm;
    Y3 = f(x3);
    return x3;

}