Метод кубической интерполяции

#include <iostream>
#include <cmath>
using namespace std;
void Coeff_11(double a, double b, double *f_x, double *c, double *d);
void kub (double a, double b, double *f_x, double *c, double *d, double *z);
double F(double x);
double F_D(double x);
 
int main () {
    double a=0,b=5;// Начало и конец отрезка
    double f_x[2], c[2], d[2], z[2];  // f_x - значение функции на границах отрезка; c - коэффициенты; d - значение производной на концах отрезка; z - экстремумы построенного многочлена
    kub (a,b,f_x,c,d,z); 
}
 
void kub (double a, double b, double *f_x, double *c, double *d, double *z){
    double min;
    f_x[0] = F(a); f_x[1] = F(b); // Заполняю массивы
    d[0] = F_D(a); d[1] = F_D(b);
    if (F_D(a)<=0 && F_D(b)>=0){ // Проверяю условие применимости метода
        Coeff_11(a,b,f_x,c,d); // Вызываю функцию вычисления 
        z[0] = (-c[2] + sqrt (c[2]*c[2] - 3*c[3]*c[1]))/(3*c[3]); // Вычисляю нули производной 
        z[1] = (-c[2] - sqrt (c[2]*c[2] - 3*c[3]*c[1]))/(3*c[3]);
        if (z[0]>=a && z[0]<=b) min = z[0]; // Если нуль входит в наш отрезок, то он является минимумом
        else min = z[1];
        if (fabs(F_D(min))<0.001) {cout<<(a + b)/2<<endl; return;} //Если значение производной исходной функции в точке найденного минимума достаточно мало, то мы нашли отрезок, в котором находится минимум исходной функции
        if (F_D(min)>0) {kub (a,min,f_x,c,d,z);} // Если значение производной исходной функции в точке найденного минимума больше нуля, то отрезком локализации минимума является отрезок (a,min)
        if (F_D(min)<0) {kub (min,b,f_x,c,d,z);} // Если значение производной исходной функции в точке найденного минимума меньше нуля, то отрезком локализации минимума является отрезок (min,b)
    }
    else cout<<"ERROR"; //Если условие применимости метода не выполняется, то вывести на экран ERROR
}
 
void Coeff_11(double a, double b, double *f_x, double *c, double *d)
{
    double tmp1, tmp2;
    //Тут считаются коэффициенты
    c[0] = f_x[0];
    c[1] = d[0];
    tmp2 = b - a;
    tmp1 = (f_x[1] - f_x[0]) / tmp2;
    c[2] = (3 * tmp1 - 2 * d[0] - d[1]) / tmp2;
    tmp2 *= tmp2;
    c[3] = (d[0] + d[1] - 2 * tmp1) / tmp2;
}
 
double F(double x) {return 3*x*x*x*x*x*x + (x-3)*(x-3);} // Функция
double F_D(double x){return 2 * (9*x*x*x*x*x + x - 3);} //Производная функции